Whitespace cleanups mostly involved removing trailing whitespace.
-#!/usr/bin/env python
-
'''
csvdriver.py
Simple driver to read general comma-separated values in Hooke
Columns are read this way:
-
+
X1 , Y1 , X2 , Y2 , X3 , Y3 ...
If the number of columns is odd, the last column is ignored.
import csv
class csvdriverDriver(lhc.Driver):
-
+
def __init__(self, filename):
-
+
self.filedata = open(filename,'r')
- self.data = list(self.filedata)
+ self.data = list(self.filedata)
self.filedata.close()
-
+
self.filetype = 'generic'
self.experiment = ''
-
+
self.filename=filename
-
+
def is_me(self):
myfile=file(self.filename)
headerline=myfile.readlines()[0]
myfile.close()
-
+
#using a custom header makes things much easier...
#(looking for raw CSV data is at strong risk of confusion)
if headerline[:-1]=='Hooke data':
return True
else:
return False
-
+
def close_all(self):
self.filedata.close()
-
+
def default_plots(self):
rrows=csv.reader(self.data)
rows=list(rrows) #transform the csv.reader iterator in a normal list
columns=lh.transposed2(rows[1:])
-
+
main_plot=lhc.PlotObject()
main_plot.vectors=[]
-
+
for index in range(0,len(columns),2):
main_plot.vectors.append([])
temp_x=columns[index]
temp_y=columns[index+1]
-
+
#convert to float (the csv gives strings)
temp_x=[float(item) for item in temp_x]
temp_y=[float(item) for item in temp_y]
-
+
main_plot.vectors[-1].append(temp_x)
main_plot.vectors[-1].append(temp_y)
-
+
main_plot.units=['x','y']
main_plot.title=self.filename
main_plot.destination=0
-
+
return [main_plot]
-
-
\ No newline at end of file
-#!/usr/bin/env python
-
'''
libhemingclamp.py
Library for interpreting Hemingway force spectroscopy files.
-Copyright (C) 2008 Massimo Sandal, Marco Brucale (University of Bologna, Italy)
+Copyright (C) 2008 Massimo Sandal, Marco Brucale (University of Bologna, Italy)
This program is released under the GNU General Public License version 2.
'''
2007_02_07: Initial implementation
'''
import string
-import libhookecurve as lhc
+import libhookecurve as lhc
class DataChunk(list):
'''Dummy class to provide ext and ret methods to the data list.
In this case ext and self can be equal.
'''
-
+
def ext(self):
return self
-
+
def ret(self):
return self
class hemingclampDriver(lhc.Driver):
-
+
def __init__(self, filename):
-
+
self.filedata = open(filename,'r')
self.data = self.filedata.readlines()[6:]
self.filedata.close()
-
+
self.filetype = 'hemingclamp'
self.experiment = 'clamp'
-
+
self.filename=filename
-
+
def __del__(self):
- self.filedata.close()
-
+ self.filedata.close()
+
def is_me(self):
'''
we define our magic heuristic for HemingClamp files
return True
else:
return False
-
+
def _getdata_all(self):
time = []
phase = []
imposed = []
trim_indexes = []
trim_counter = 0.0
-
+
for i in self.data:
temp = string.split(i)
#time.append(float(temp[0])*(1.0e-3)) # This is managed differently now, since each data point = 1ms: see below
if phase[x] != trim_counter:
trim_indexes.append(x)
trim_counter = phase[x]
-
+
#we rebuild the time counter assuming 1 point = 1 millisecond
c=0.0
for z in zpiezo:
time.append(c)
- c+=(1.0e-3)
-
+ c+=(1.0e-3)
+
return time,phase,zpiezo,defl,imposed,trim_indexes
-
+
def time(self):
return DataChunk(self._getdata_all()[0])
def phase(self):
return DataChunk(self._getdata_all()[1])
-
+
def zpiezo(self):
return DataChunk(self._getdata_all()[2])
-
+
def deflection(self):
return DataChunk(self._getdata_all()[3])
def trimindexes(self):
return DataChunk(self._getdata_all()[5])
-
+
def close_all(self):
'''
Explicitly closes all files
'''
self.filedata.close()
-
+
def default_plots(self):
main_plot=lhc.PlotObject()
defl_plot=lhc.PlotObject()
-
+
time=self.time()
phase=self.phase()
zpiezo=self.zpiezo()
deflection=self.deflection()
imposed=self.imposed()
-
+
main_plot.vectors=[[time,zpiezo],[time,phase]]
main_plot.units=['seconds','meters']
main_plot.destination=0
main_plot.title=self.filename
-
+
defl_plot.vectors=[[time,deflection],[time,imposed]]
defl_plot.units=['seconds','Newtons']
defl_plot.destination=1
-
+
return [main_plot, defl_plot]
-
\ No newline at end of file
+
-#!/usr/bin/env python
-
import string
-import libhookecurve as lhc
+import libhookecurve as lhc
class DataChunk(list):
#Dummy class to provide ext and ret methods to the data list.
-
+
def ext(self):
halflen=(len(self)/2)
return self[0:halflen]
-
+
def ret(self):
halflen=(len(self)/2)
return self[halflen:]
self.filename=filename #self.filename can always be useful, and should be defined
self.filedata = open(filename,'r') #We open the file
self.filelines=self.filedata.readlines()
- self.filedata.close()
+ self.filedata.close()
'''These are two strings that can be used by Hooke commands/plugins to understand what they are looking at. They have no other
meaning. They have to be somehow defined however - commands often look for those variables.
-
+
self.filetype should contain the name of the exact filetype defined by the driver (so that filetype-specific commands can know
if they're dealing with the correct filetype)
self.experiment should contain instead the type of data involved (for example, various drivers can be used for force-clamp experiments,
- but hooke commands could like to know if we're looking at force clamp data, regardless of their origin, and not other
+ but hooke commands could like to know if we're looking at force clamp data, regardless of their origin, and not other
kinds of data)
-
+
Of course, all other variables you like can be defined in the class.
'''
self.filetype = 'jpk'
self.experiment = 'smfs'
-
-
+
+
def __del__(self):
- self.filedata.close()
-
+ self.filedata.close()
+
def is_me(self):
'''
we define our magic heuristic for jpk files
return True
else:
return False
-
+
def close_all(self):
self.filedata.close()
-
+
def _read_data_segment(self):
#routine that actually reads the data
-
+
height_ms=[]
height_m=[]
height=[]
v_deflection=[]
h_deflection=[]
-
+
self.springconstant=0 #if we don't meet any spring constant, use deflection...
-
+
for line in self.filelines:
#we meet the segment defining the order of data columns
-
+
if line[0:9]=='# columns':
splitline=line.split()[2:]
height_ms_index=splitline.index('smoothedStrainGaugeHeight')
height_index=splitline.index('height')
v_deflection_index=splitline.index('vDeflection')
#h_deflection=splitline.index('hDeflection')
-
+
if line[0:16]=='# springConstant':
self.springconstant=float(line.split()[2])
-
+
if line[0] != '#' and len(line.split())>1:
dataline=line.split()
height_ms.append(float(dataline[height_ms_index]))
height.append(float(dataline[height_index]))
v_deflection.append(float(dataline[v_deflection_index]))
#h_deflection.append(float(dataline[h_deflection_index]))
-
+
if self.springconstant != 0:
force=[item*self.springconstant for item in v_deflection]
else: #we have measured no spring constant :(
force=v_deflection
-
+
height_ms=DataChunk([item*-1 for item in height_ms])
height_m=DataChunk([item*-1 for item in height_m])
height=DataChunk([item*-1 for item in height])
deflection=DataChunk(v_deflection)
force=DataChunk(force)
-
+
return height_ms,height_m,height,deflection,force
-
+
def deflection(self):
height_ms,height_m,height,deflection,force=self._read_data_segment()
deflection_ext=deflection.ext()
deflection_ret=deflection.ret()
deflection_ret.reverse()
return deflection_ext,deflection_ret
-
+
def default_plots(self):
-
+
height_ms,height_m,height,deflection,force=self._read_data_segment()
-
+
height_ms_ext=height_ms.ext()
height_ms_ret=height_ms.ret()
force_ext=force.ext()
#reverse the return data, to make it coherent with hooke standard
height_ms_ret.reverse()
force_ret.reverse()
-
- main_plot=lhc.PlotObject()
+
+ main_plot=lhc.PlotObject()
main_plot.add_set(height_ms_ext,force_ext)
main_plot.add_set(height_ms_ret,force_ret)
-
-
-
+
+
+
if self.springconstant != 0:
main_plot.units=['meters','force']
else:
main_plot.units=['meters','meters']
-
+
main_plot.normalize_vectors()
-
+
main_plot.destination=0
main_plot.title=self.filename
-
- return [main_plot]
\ No newline at end of file
+
+ return [main_plot]
-#!/usr/bin/env python
-
'''
mcs.py
import struct
class mcsDriver(lhc.Driver):
-
+
def __init__(self, filename):
'''
Open the RED (A) ones; the BLUE (D) mirror ones will be automatically opened
oth[-8]='d'
othername=''.join(oth)
self.filename=filename
- self.othername=othername
-
+ self.othername=othername
+
#print self.filename, self.othername
-
+
self.filedata=open(filename,'rb')
self.reddata=self.filedata.read()
self.filedata.close()
-
+
self.filebluedata=open(othername,'rb') #open also the blue ones
self.bluedata=self.filebluedata.read()
self.filebluedata.close()
-
+
self.filetype = 'mcs'
self.experiment = 'smfluo'
-
+
def is_me(self):
if self.filename[-3:].lower()=='mcs':
return True
else:
return False
-
+
def close_all(self):
self.filedata.close()
self.filebluedata.close()
-
-
+
+
def default_plots(self):
red_data=self.read_file(self.reddata)
blue_data=self.read_file(self.bluedata)
blue_data=[-1*float(item) for item in blue_data] #visualize blue as "mirror" of red
-
+
main_plot=lhc.PlotObject()
main_plot.add_set(range(len(red_data)),red_data)
main_plot.add_set(range(len(blue_data)),blue_data)
main_plot.destination=0
main_plot.title=self.filename
main_plot.colors=['red','blue']
-
+
return [main_plot]
-
- def read_file(self, raw_data):
+
+ def read_file(self, raw_data):
real_data=[]
intervalsperfile=struct.unpack('h', raw_data[10:12])[0] #read in number of intervals in this file
#this data is contained in bit offset 10-12 in mcs file
#see http://docs.python.org/library/struct.html#module-struct for additional explanation
-
+
numbytes=len(raw_data) #data is stored in 4-byte chunks, starting with pos 256
for j in range(0,intervalsperfile): #read in all intervals in file
temp=raw_data[256+j*4:256+j*4+4] #data starts at byte offset 256
real_data.append(struct.unpack('i', temp)[0]) #[0] because it returns a 1-element tuple
- return real_data
\ No newline at end of file
+ return real_data
-#!/usr/bin/env python
-
'''
mfp1dexport.py
import libhooke as lh
class mfp1dexportDriver(lhc.Driver):
-
+
def __init__(self, filename):
-
+
self.filename=filename
self.filedata=open(filename,'rU')
self.lines=list(self.filedata.readlines())
self.filedata.close()
-
+
self.filetype='mfp1dexport'
self.experiment='smfs'
-
+
def close_all(self):
self.filedata.close()
-
+
def is_me(self):
try:
self.raw_header=self.lines[0:38]
except:
#Not enough lines for a header; not a good file
return False
-
+
#FIXME: We want a more reasonable header recognition
if self.raw_header[0][0:4]=='Wave':
return True
else:
return False
-
+
def _read_columns(self):
-
+
self.raw_columns=self.lines[39:]
-
+
kline=None
for line in self.lines:
if line[:7]=='SpringC':
kline=line
break
-
+
kline=kline.split(':')
-
+
#self.k=float(self.raw_header[23][8:])
self.k=float(kline[1])
-
-
+
+
xext=[]
xret=[]
yext=[]
yext.append(float(spline[1]))
xret.append(float(spline[2]))
yret.append(float(spline[3]))
-
+
return [[xext,yext],[xret,yret]]
-
+
def deflection(self):
self.data=self._read_columns()
return self.data[0][1],self.data[1][1]
-
-
- def default_plots(self):
+
+
+ def default_plots(self):
main_plot=lhc.PlotObject()
defl_ext,defl_ret=self.deflection()
yextforce=[i*self.k for i in defl_ext]
-#!/usr/bin/env python
-
'''
libpicoforce.py
class DataChunk(list):
#Dummy class to provide ext and ret methods to the data list.
-
+
def ext(self):
halflen=(len(self)/2)
return self[0:halflen]
-
+
def ret(self):
halflen=(len(self)/2)
return self[halflen:]
class picoforceDriver(lhc.Driver):
#Construction and other special methods
-
+
def __init__(self,filename):
'''
constructor method
'''
-
+
self.textfile=file(filename)
self.binfile=file(filename,'rb')
-
+
#The 0,1,2 data chunks are:
#0: D (vs T)
#1: Z (vs T)
#2: D (vs Z)
-
-
+
+
self.filepath=filename
self.debug=False
-
+
self.filetype='picoforce'
self.experiment='smfs'
-
-
+
+
#Hidden methods. These are meant to be used only by API functions. If needed, however,
#they can be called just like API methods.
-
+
def _get_samples_line(self):
'''
Gets the samples per line parameters in the file, to understand trigger behaviour.
'''
self.textfile.seek(0)
-
+
samps_expr=re.compile(".*Samps")
-
+
samps_values=[]
for line in self.textfile.readlines():
if samps_expr.match(line):
samps_values.append(samps)
except:
pass
-
+
#We raise a flag for the fact we meet an offset, otherwise we would take spurious data length arguments.
-
+
return int(samps_values[0])
-
+
def _get_chunk_coordinates(self):
'''
This method gets the coordinates (offset and length) of a data chunk in our
Picoforce file.
-
- It returns a list containing two tuples:
- the first element of each tuple is the data_offset, the second is the corresponding
+
+ It returns a list containing two tuples:
+ the first element of each tuple is the data_offset, the second is the corresponding
data size.
-
- In near future probably each chunk will get its own data structure, with
+
+ In near future probably each chunk will get its own data structure, with
offset, size, type, etc.
'''
self.textfile.seek(0)
-
+
offset_expr=re.compile(".*Data offset")
length_expr=re.compile(".*Data length")
offset=int(line.split()[2]) #the third word splitted is the offset (in bytes)
data_offsets.append(offset)
#We raise a flag for the fact we meet an offset, otherwise we would take spurious data length arguments.
- flag_offset=1
-
+ flag_offset=1
+
#same for the data length
- if length_expr.match(line) and flag_offset:
+ if length_expr.match(line) and flag_offset:
size=int(line.split()[2])
data_sizes.append(size)
#Put down the offset flag until the next offset is met.
flag_offset=0
return zip(data_offsets,data_sizes)
-
+
def _get_data_chunk(self,whichchunk):
'''
reads a data chunk and converts it in 16bit signed int.
'''
offset,size=self._get_chunk_coordinates()[whichchunk]
-
-
+
+
self.binfile.seek(offset)
raw_chunk=self.binfile.read(size)
-
+
my_chunk=[]
for data_position in range(0,len(raw_chunk),2):
data_unit_bytes=raw_chunk[data_position:data_position+2]
#The unpack function converts 2-bytes in a signed int ('h').
#we use output[0] because unpack returns a 1-value tuple, and we want the number only
data_unit=struct.unpack('h',data_unit_bytes)[0]
- my_chunk.append(data_unit)
-
+ my_chunk.append(data_unit)
+
return DataChunk(my_chunk)
-
+
def _get_Zscan_info(self,index):
'''
gets the Z scan informations needed to interpret the data chunk.
These info come from the general section, BEFORE individual chunk headers.
-
+
By itself, the function will parse for three parameters.
(index) that tells the function what to return when called by
exposed API methods.
index=2 : returns Zscan_V_size
'''
self.textfile.seek(0)
-
+
ciaoforcelist_expr=re.compile(".*Ciao force")
zscanstart_expr=re.compile(".*@Z scan start")
zscansize_expr=re.compile(".*@Z scan size")
-
+
ciaoforce_flag=0
theline=0
for line in self.textfile.readlines():
if ciaoforcelist_expr.match(line):
ciaoforce_flag=1 #raise a flag: zscanstart and zscansize params to read are later
-
+
if ciaoforce_flag and zscanstart_expr.match(line):
raw_Zscanstart_line=line.split()
-
+
if ciaoforce_flag and zscansize_expr.match(line):
raw_Zscansize_line=line.split()
for itemscanstart,itemscansize in zip(raw_Zscanstart_line,raw_Zscansize_line):
Zscanstart_line.append(itemscanstart.strip('[]()'))
Zscansize_line.append(itemscansize.strip('[]()'))
-
+
Zscan_V_LSB=float(Zscanstart_line[6])
Zscan_V_start=float(Zscanstart_line[8])
Zscan_V_size=float(Zscansize_line[8])
-
+
return (Zscan_V_LSB,Zscan_V_start,Zscan_V_size)[index]
-
+
def _get_Z_magnify_scale(self,whichchunk):
'''
gets Z scale and Z magnify
Here we get Z scale/magnify from the 'whichchunk' only.
whichchunk=1,2,3
TODO: make it coherent with data_chunks syntaxis (0,1,2)
-
+
In future, should we divide the *file* itself into chunk descriptions and gain
true chunk data structures?
'''
self.textfile.seek(0)
-
+
z_scale_expr=re.compile(".*@4:Z scale")
z_magnify_expr=re.compile(".*@Z magnify")
-
+
ramp_size_expr=re.compile(".*@4:Ramp size")
ramp_offset_expr=re.compile(".*@4:Ramp offset")
-
+
occurrences=0
found_right=0
-
-
+
+
for line in self.textfile.readlines():
if z_magnify_expr.match(line):
occurrences+=1
raw_z_magnify_expression=line.split()
else:
found_right=0
-
+
if found_right and z_scale_expr.match(line):
raw_z_scale_expression=line.split()
if found_right and ramp_size_expr.match(line):
raw_ramp_size_expression=line.split()
if found_right and ramp_offset_expr.match(line):
raw_ramp_offset_expression=line.split()
-
+
return float(raw_z_magnify_expression[5]),float(raw_z_scale_expression[7]), float(raw_ramp_size_expression[7]), float(raw_ramp_offset_expression[7]), float(raw_z_scale_expression[5][1:])
-
-
+
+
#Exposed APIs.
- #These are the methods that are meant to be called from external apps.
-
+ #These are the methods that are meant to be called from external apps.
+
def LSB_to_volt(self,chunknum,voltrange=20):
'''
Converts the LSB data of a given chunk (chunknum=0,1,2) in volts.
First step to get the deflection and the force.
-
+
SYNTAXIS:
item.LSB_to_volt(chunknum, [voltrange])
-
+
The voltrange is by default set to 20 V.
'''
return DataChunk([((float(lsb)/65535)*voltrange) for lsb in self.data_chunks[chunknum]])
-
+
def LSB_to_deflection(self,chunknum,deflsensitivity=None,voltrange=20):
'''
Converts the LSB data in deflection (meters).
-
+
SYNTAXIS:
item.LSB_to_deflection(chunknum, [deflection sensitivity], [voltrange])
-
+
chunknum is the chunk you want to parse (0,1,2)
-
- The deflection sensitivity by default is the one parsed from the file.
+
+ The deflection sensitivity by default is the one parsed from the file.
The voltrange is by default set to 20 V.
'''
if deflsensitivity is None:
deflsensitivity=self.get_deflection_sensitivity()
-
- lsbvolt=self.LSB_to_volt(chunknum)
+
+ lsbvolt=self.LSB_to_volt(chunknum)
return DataChunk([volt*deflsensitivity for volt in lsbvolt])
-
+
def deflection(self):
'''
Get the actual force curve deflection.
'''
deflchunk= self.LSB_to_deflection(2)
return deflchunk.ext(),deflchunk.ret()
-
+
def LSB_to_force(self,chunknum=2,Kspring=None,voltrange=20):
'''
Converts the LSB data (of deflection) in force (newtons).
-
+
SYNTAXIS:
item.LSB_to_force([chunknum], [spring constant], [voltrange])
-
+
chunknum is the chunk you want to parse (0,1,2). The chunk used is by default 2.
The spring constant by default is the one parsed from the file.
The voltrange is by default set to 20 V.
'''
if Kspring is None:
Kspring=self.get_spring_constant()
-
- lsbdefl=self.LSB_to_deflection(chunknum)
+
+ lsbdefl=self.LSB_to_deflection(chunknum)
return DataChunk([(meter*Kspring) for meter in lsbdefl])
-
+
def get_Zscan_V_start(self):
return self._get_Zscan_info(1)
-
+
def get_Zscan_V_size(self):
return self._get_Zscan_info(2)
-
+
def get_Z_scan_sensitivity(self):
'''
gets Z sensitivity
'''
self.textfile.seek(0)
-
+
z_sensitivity_expr=re.compile(".*@Sens. Zsens")
-
+
for line in self.textfile.readlines():
if z_sensitivity_expr.match(line):
z_sensitivity=float(line.split()[3])
#return it in SI units (that is: m/V, not nm/V)
return z_sensitivity*(10**(-9))
-
+
def get_Z_magnify(self,whichchunk):
'''
Gets the Z magnify factor. Normally it is 1, unknown exact use as of 2006-01-13
'''
return self._get_Z_magnify_scale(whichchunk)[0]
-
+
def get_Z_scale(self,whichchunk):
'''
Gets the Z scale.
'''
return self._get_Z_magnify_scale(whichchunk)[1]
-
+
def get_ramp_size(self,whichchunk):
'''
Gets the -user defined- ramp size
'''
return self._get_Z_magnify_scale(whichchunk)[2]
-
+
def get_ramp_offset(self,whichchunk):
'''
Gets the ramp offset
'''
return self._get_Z_magnify_scale(whichchunk)[3]
-
+
def get_Z_scale_LSB(self,whichchunk):
'''
Gets the LSB-to-volt conversion factor of the Z data.
(so called hard-scale in the Nanoscope documentation)
-
+
'''
return self._get_Z_magnify_scale(whichchunk)[4]
-
+
def get_deflection_sensitivity(self):
'''
gets deflection sensitivity
- '''
+ '''
self.textfile.seek(0)
-
+
def_sensitivity_expr=re.compile(".*@Sens. DeflSens")
-
+
for line in self.textfile.readlines():
if def_sensitivity_expr.match(line):
def_sensitivity=float(line.split()[3])
break
#return it in SI units (that is: m/V, not nm/V)
return def_sensitivity*(10**(-9))
-
+
def get_spring_constant(self):
'''
gets spring constant.
They are normally all equal, but we retain all three for future...
'''
self.textfile.seek(0)
-
+
springconstant_expr=re.compile(".*Spring Constant")
-
+
constants=[]
-
+
for line in self.textfile.readlines():
if springconstant_expr.match(line):
constants.append(float(line.split()[2]))
-
+
return constants[0]
-
+
def get_Zsensorsens(self):
'''
gets Zsensorsens for Z data.
-
+
This is the sensitivity needed to convert the LSB data in nanometers for the Z-vs-T data chunk.
- '''
+ '''
self.textfile.seek(0)
-
+
zsensorsens_expr=re.compile(".*Sens. ZSensorSens")
-
+
for line in self.textfile.readlines():
if zsensorsens_expr.match(line):
zsensorsens_raw_expression=line.split()
#we must take only first occurrence, so we exit from the cycle immediately
break
-
+
return (float(zsensorsens_raw_expression[3]))*(10**(-9))
-
+
def Z_data(self):
'''
returns converted ext and ret Z curves.
#rampsize_zt=self.get_ramp_size(2)
#rampoffset_zt=self.get_ramp_offset(2)
zsensorsens=self.get_Zsensorsens()
-
+
'''
The magic formula that converts the Z data is:
-
+
meters = LSB * V_lsb_conversion_factor * ZSensorSens
'''
-
+
#z_curves=[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].pair['ext']],[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].pair['ret']]
- z_curves=[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].ext()],[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].ret()]
+ z_curves=[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].ext()],[item*Zlsb_zt*zsensorsens for item in self.data_chunks[1].ret()]
z_curves=[DataChunk(item) for item in z_curves]
return z_curves
-
+
def Z_extremes(self):
'''
returns the extremes of the Z values
z_extremes={}
z_extremes['ext']=zcurves[0][0],zcurves[0][-1]
z_extremes['ret']=zcurves[1][0],zcurves[1][-1]
-
+
return z_extremes
-
+
def Z_step(self):
'''
returns the calculated step between the Z values
zpoints={}
z_extremes=self.Z_extremes()
-
+
zrange['ext']=abs(z_extremes['ext'][0]-z_extremes['ext'][1])
zrange['ret']=abs(z_extremes['ret'][0]-z_extremes['ret'][1])
-
+
#We must take 1 from the calculated zpoints, or when I use the arange function gives me a point more
#with the step. That is, if I have 1000 points, and I use arange(start,stop,step), I have 1001 points...
#For cleanness, solution should really be when using arange, but oh well...
zpoints['ret']=len(self.Z_data()[1])-1
#this syntax must become coherent!!
return (zrange['ext']/zpoints['ext']),(zrange['ret']/zpoints['ret'])
-
+
def Z_domains(self):
'''
returns the Z domains on which to plot the force data.
-
+
The Z domains are returned as a single long DataChunk() extended list. The extension and retraction part
- can be extracted using ext() and ret() methods.
- '''
+ can be extracted using ext() and ret() methods.
+ '''
x1step=self.Z_step()[0]
- x2step=self.Z_step()[1]
-
+ x2step=self.Z_step()[1]
+
try:
xext=arange(self.Z_extremes()['ext'][0],self.Z_extremes()['ext'][1],-x1step)
xret=arange(self.Z_extremes()['ret'][0],self.Z_extremes()['ret'][1],-x2step)
xext=arange(0,1)
xret=arange(0,1)
print 'picoforce.py: Warning. xext, xret domains cannot be extracted.'
-
+
if not (len(xext)==len(xret)):
if self.debug:
#print warning
print "You cannot trust the resulting curve."
print "Until a solution is found, I substitute the ext domain with the ret domain. Sorry."
xext=xret
-
+
return DataChunk(xext.tolist()+xret.tolist())
-
+
def Z_scan_size(self):
return self.get_Zscan_V_size()*self.get_Z_scan_sensitivity()
-
+
def Z_start(self):
return self.get_Zscan_V_start()*self.get_Z_scan_sensitivity()
-
+
def ramp_size(self,whichchunk):
'''
to be implemented if needed
'''
raise "Not implemented yet."
-
-
+
+
def ramp_offset(self,whichchunk):
'''
to be implemented if needed
'''
raise "Not implemented yet."
-
+
def detriggerize(self, forcext):
'''
Cuts away the trigger-induced s**t on the extension curve.
DEPRECATED
cutindex=2
startvalue=forcext[0]
-
- for index in range(len(forcext)-1,2,-2):
+
+ for index in range(len(forcext)-1,2,-2):
if forcext[index]>startvalue:
cutindex=index
else:
return cutindex
'''
return 0
-
+
def is_me(self):
'''
self-identification of file type magic
curve_file=file(self.filepath)
header=curve_file.read(30)
curve_file.close()
-
+
if header[2:17] == 'Force file list': #header of a picoforce file
self.data_chunks=[self._get_data_chunk(num) for num in [0,1,2]]
return True
else:
return False
-
+
def close_all(self):
'''
Explicitly closes all files
'''
self.textfile.close()
self.binfile.close()
-
+
def default_plots(self):
'''
creates the default PlotObject
'''
-
-
+
+
force=self.LSB_to_force()
zdomain=self.Z_domains()
-
+
samples=self._get_samples_line()
#cutindex=0
#cutindex=self.detriggerize(force.ext())
-
+
main_plot=lhc.PlotObject()
-
+
main_plot.vectors=[[zdomain.ext()[0:samples], force.ext()[0:samples]],[zdomain.ret()[0:samples], force.ret()[0:samples]]]
main_plot.normalize_vectors()
main_plot.units=['meters','newton']
main_plot.destination=0
main_plot.title=self.filepath
-
-
+
+
return [main_plot]
-#!/usr/bin/env python
-
'''
libpicoforce.py
class DataChunk(list):
#Dummy class to provide ext and ret methods to the data list.
-
+
def ext(self):
halflen=(len(self)/2)
return self[0:halflen]
-
+
def ret(self):
halflen=(len(self)/2)
return self[halflen:]
class picoforcealtDriver(lhc.Driver):
#Construction and other special methods
-
+
def __init__(self,filename):
'''
constructor method
'''
-
+
self.textfile=file(filename)
self.binfile=file(filename,'rb')
-
+
#The 0,1,2 data chunks are:
#0: D (vs T)
#1: Z (vs T)
#2: D (vs Z)
-
+
self.forcechunk=0
self.distancechunk=1
#TODO eliminate the need to set chunk numbers
-
+
self.filepath=filename
self.debug=True
-
+
self.filetype='picoforce'
self.experiment='smfs'
-
-
-
+
+
+
def _get_samples_line(self):
'''
Gets the samples per line parameters in the file, to understand trigger behaviour.
'''
self.textfile.seek(0)
-
+
samps_expr=re.compile(".*Samps")
-
+
samps_values=[]
for line in self.textfile.readlines():
if samps_expr.match(line):
samps_values.append(samps)
except:
pass
-
+
#We raise a flag for the fact we meet an offset, otherwise we would take spurious data length arguments.
-
+
return int(samps_values[0])
-
+
def _get_chunk_coordinates(self):
'''
This method gets the coordinates (offset and length) of a data chunk in our
Picoforce file.
-
- It returns a list containing two tuples:
- the first element of each tuple is the data_offset, the second is the corresponding
+
+ It returns a list containing two tuples:
+ the first element of each tuple is the data_offset, the second is the corresponding
data size.
-
- In near future probably each chunk will get its own data structure, with
+
+ In near future probably each chunk will get its own data structure, with
offset, size, type, etc.
'''
self.textfile.seek(0)
-
+
offset_expr=re.compile(".*Data offset")
length_expr=re.compile(".*Data length")
offset=int(line.split()[2]) #the third word splitted is the offset (in bytes)
data_offsets.append(offset)
#We raise a flag for the fact we meet an offset, otherwise we would take spurious data length arguments.
- flag_offset=1
-
+ flag_offset=1
+
#same for the data length
- if length_expr.match(line) and flag_offset:
+ if length_expr.match(line) and flag_offset:
size=int(line.split()[2])
data_sizes.append(size)
#Put down the offset flag until the next offset is met.
flag_offset=0
return zip(data_offsets,data_sizes)
-
+
def _get_data_chunk(self,whichchunk):
'''
reads a data chunk and converts it in 16bit signed int.
'''
offset,size=self._get_chunk_coordinates()[whichchunk]
-
-
+
+
self.binfile.seek(offset)
raw_chunk=self.binfile.read(size)
-
+
my_chunk=[]
for data_position in range(0,len(raw_chunk),2):
data_unit_bytes=raw_chunk[data_position:data_position+2]
#The unpack function converts 2-bytes in a signed int ('h').
#we use output[0] because unpack returns a 1-value tuple, and we want the number only
data_unit=struct.unpack('h',data_unit_bytes)[0]
- my_chunk.append(data_unit)
-
+ my_chunk.append(data_unit)
+
return DataChunk(my_chunk)
def _force(self):
voltrange=1
z_scale=self._get_Z_scale()
deflsensitivity=self.get_deflection_sensitivity()
- volts=[((float(lsb))*voltrange*z_scale) for lsb in self.data_chunks[self.forcechunk]]
+ volts=[((float(lsb))*voltrange*z_scale) for lsb in self.data_chunks[self.forcechunk]]
deflect=[volt*deflsensitivity for volt in volts]
-
+
return deflect
-
-
+
+
def _Z(self):
- #returns distance vector (calculated instead than from data chunk)
+ #returns distance vector (calculated instead than from data chunk)
rampsize=self._get_rampsize()
sampsline=self._get_samples_line()
senszscan=self._get_Z_scan_sens()
-
+
xstep=senszscan*rampsize/sampsline*10**(-9)
xext=arange(sampsline*xstep,0,-xstep)
xret=arange(sampsline*xstep,0,-xstep)
-
+
return DataChunk(xext.tolist()+xret.tolist())
-
+
def _get_Z_scale(self):
self.textfile.seek(0)
expr=re.compile(".*@4:Z scale")
-
+
for line in self.textfile.readlines():
if expr.match(line):
zscale=float((line.split()[5]).strip("() []"))
break
return zscale
-
+
def _get_rampsize(self):
self.textfile.seek(0)
expr=re.compile(".*@4:Ramp size:")
-
+
for line in self.textfile.readlines():
if expr.match(line):
zsens=float((line.split()[7]).strip("() []"))
break
return zsens
-
+
def _get_Z_scan_sens(self):
self.textfile.seek(0)
expr=re.compile(".*@Sens. Zsens")
-
+
for line in self.textfile.readlines():
if expr.match(line):
zsens=float((line.split()[3]).strip("() []"))
break
return zsens
-
-
-
+
+
+
def get_deflection_sensitivity(self):
'''
gets deflection sensitivity
- '''
+ '''
self.textfile.seek(0)
-
+
def_sensitivity_expr=re.compile(".*@Sens. DeflSens")
-
+
for line in self.textfile.readlines():
if def_sensitivity_expr.match(line):
def_sensitivity=float(line.split()[3])
break
#return it in SI units (that is: m/V, not nm/V)
return def_sensitivity*(10**(-9))
-
+
def get_spring_constant(self):
'''
gets spring constant.
They are normally all equal, but we retain all three for future...
'''
self.textfile.seek(0)
-
+
springconstant_expr=re.compile(".*Spring Constant")
-
+
constants=[]
-
+
for line in self.textfile.readlines():
if springconstant_expr.match(line):
constants.append(float(line.split()[2]))
-
+
return constants[0]
-
+
def is_me(self):
'''
self-identification of file type magic
curve_file=file(self.filepath)
header=curve_file.read(30)
curve_file.close()
-
+
if header[2:17] == 'Force file list': #header of a picoforce file
#here DONT translate chunk
self.data_chunks=[self._get_data_chunk(num) for num in [0,1,2]]
return True
else:
return False
-
+
def close_all(self):
'''
Explicitly closes all files
'''
self.textfile.close()
self.binfile.close()
-
+
def default_plots(self):
'''
creates the default PlotObject
main_plot.units=['meters','newton']
main_plot.destination=0
main_plot.title=self.filepath
-
-
+
+
return [main_plot]
def deflection(self):
-#!/usr/bin/env python
-
'''
tutorialdriver.py
'''
This is a *generic* driver, not a specific force spectroscopy driver.
See the written documentation to see what a force spectroscopy driver must be defined to take into account Hooke facilities.
-
+
Our driver is a class with the name convention nameofthedriverDriver, where "nameofthedriver" is the filename.
The driver must inherit from the parent class lhc.Driver, so the syntax is
class nameofthedriverDriver(lhc.Driver)
'''
In this case, we have a data format that is just a list of ASCII values, so we can just divide that in rows, and generate a list
with each item being a row.
- Of course if your data files are binary, or follow a different approach, do whatever you need. :)
+ Of course if your data files are binary, or follow a different approach, do whatever you need. :)
'''
self.data = list(self.filedata)
self.filedata.close() #remember to close the file
-
+
'''These are two strings that can be used by Hooke commands/plugins to understand what they are looking at. They have no other
meaning. They have to be somehow defined however - commands often look for those variables.
-
+
self.filetype should contain the name of the exact filetype defined by the driver (so that filetype-specific commands can know
if they're dealing with the correct filetype)
self.experiment should contain instead the type of data involved (for example, various drivers can be used for force-clamp experiments,
- but hooke commands could like to know if we're looking at force clamp data, regardless of their origin, and not other
+ but hooke commands could like to know if we're looking at force clamp data, regardless of their origin, and not other
kinds of data)
-
+
Of course, all other variables you like can be defined in the class.
'''
self.filetype = 'tutorial'
self.experiment = 'generic'
-
+
def is_me(self):
'''
THIS METHOD MUST BE DEFINED.
RETURNS: Boolean (True or False)
- This method must be an heuristic that looks at the file content and decides if the file can be opened by the driver itself.
- It returns True if the file opened can be interpreted by the current driver, False otherwise.
+ This method must be an heuristic that looks at the file content and decides if the file can be opened by the driver itself.
+ It returns True if the file opened can be interpreted by the current driver, False otherwise.
Defining this method allows Hooke to understand what kind of files we're looking at automatically.
-
+
We have to re-open/re-close the file here.
'''
-
+
myfile=open(self.filename, 'r')
headerline=myfile.readlines()[0] #we take the first line
myfile.close()
-
+
'''
Here, our "magic fingerprint" is the TUTORIAL_FILE header. Of course, depending on the data file, you can have interesting
headers, or patterns, etc. that you can use to guess the data format. What matters is successful recognizing, and returning
return True
else:
return False
-
+
def _generate_vectors(self):
'''
Here we parse the data and generate the raw vectors. This method has only to do with the peculiar file format here, so it's of
- no big interest (I just wrote it to present a functional driver).
-
+ no big interest (I just wrote it to present a functional driver).
+
Only thing to remember, it can be nice to call methods that are used only "internally" by the driver (or by plugins) with a
"_" prefix, so to have a visual remark. But it's just an advice.
'''
elif item[:-1]=='PLOT2':
whatplot=item[:-1]
elif item[0]=='X' or item[0]=='Y':
- pos=positions[item[:-1]]
+ pos=positions[item[:-1]]
else:
pass
-
+
return vectors
def close_all(self):
In this method, all file objects defined in the driver must be closed.
'''
self.filename.close()
-
-
+
+
def default_plots(self):
'''
THIS METHOD MUST BE DEFINED.
RETURNS: [ lhc.PlotObject ] or [ lhc.PlotObject, lhc.PlotObject]
-
+
This is the method that returns the plots to Hooke.
It must return a list with *one* or *two* PlotObjects.
-
+
See the libhookecurve.py source code to see how PlotObjects are defined and work in detail.
'''
gen_vectors=self._generate_vectors()
-
+
#Here we create the main plot PlotObject and initialize its vectors.
main_plot=lhc.PlotObject()
main_plot.vectors=[]
#The same for the other plot.
other_plot=lhc.PlotObject()
other_plot.vectors=[]
-
+
'''
Now we fill the plot vectors with our data.
set 1 set 2
The "correct" shape of the vector is [ [[x1,x2,x3...],[y1,y2,y3...]] , [[x1,x2,x3...],[y1,y2,y3...]] ], so we have to put stuff in this way into it.
-
+
The add_set() method takes care of this , just use plot.add_set(x,y).
'''
main_plot.add_set(gen_vectors['PLOT1'][0],gen_vectors['PLOT1'][1])
main_plot.add_set(gen_vectors['PLOT1'][2],gen_vectors['PLOT1'][3])
-
+
other_plot.add_set(gen_vectors['PLOT2'][0],gen_vectors['PLOT2'][1])
other_plot.add_set(gen_vectors['PLOT2'][2],gen_vectors['PLOT2'][3])
-
+
'''
normalize_vectors() trims the vectors, so that if two x/y couples are of different lengths, the latest
points are trimmed (otherwise we have a python error). Always a good idea to run it, to avoid crashes.
'''
main_plot.normalize_vectors()
other_plot.normalize_vectors()
-
+
'''
Here we define:
- units: [string, string], define the measure units of X and Y axes
- destination: 0/1 , defines where to plot the plot (0=top, 1=bottom), default=0
- title: string , the plot title.
-
+
for each plot.
Again, see libhookecurve.py comments for details.
'''
main_plot.units=['unit of x','unit of y']
main_plot.destination=0
main_plot.title=self.filename+' main'
-
+
other_plot.units=['unit of x','unit of y']
other_plot.destination=1
other_plot.title=self.filename+' other'
-
- return [main_plot, other_plot]
\ No newline at end of file
+
+ return [main_plot, other_plot]
\r
app.MainLoop()\r
\r
-main()\r
+if __name__ == '__main__':\r
+ main()\r
This program is released under the GNU General Public License version 2.
'''
-
from libhooke import * #FIXME
import libhookecurve as lhc
class HookeCli(cmd.Cmd):
-
+
def __init__(self,frame,list_of_events,events_from_gui,config,drivers):
cmd.Cmd.__init__(self)
-
+
self.prompt = 'hooke: '
-
-
+
+
self.current_list=[] #the playlist we're using
-
- self.current=None #the current curve under analysis.
+
+ self.current=None #the current curve under analysis.
self.plots=None
'''
The actual hierarchy of the "current curve" is a bit complex:
-
+
self.current = the lhc.HookeCurve container object of the current curve
self.current.curve = the current "real" curve object as defined in the filetype driver class
self.current.curve.default_plots() = the default plots of the filetype driver.
-
- The plot objects obtained by mean of self.current.curve.default_plots()
+
+ The plot objects obtained by mean of self.current.curve.default_plots()
then undergoes modifications by the plotmanip
- modifier functions. The modified plot is saved in self.plots and used if needed by other functions.
+ modifier functions. The modified plot is saved in self.plots and used if needed by other functions.
'''
-
-
+
+
self.pointer=0 #a pointer to navigate the current list
-
+
#Things that come from outside
self.frame=frame #the wx frame we refer to
self.list_of_events=list_of_events #a list of wx events we use to interact with the GUI
self.events_from_gui=events_from_gui #the Queue object we use to have messages from the GUI
self.config=config #the configuration dictionary
self.drivers=drivers #the file format drivers
-
+
#get plot manipulation functions
plotmanip_functions=[]
for object_name in dir(self):
self.plotmanip[nameindex] = item
else:
pass
-
-
+
+
self.playlist_saved=0 #Did we save the playlist?
self.playlist_name='' #Name of playlist
self.notes_saved=1 #Did we save the notes?
#create outlet
self.outlet=lout.Outlet()
-
+
#Data that must be saved in the playlist, related to the whole playlist (not individual curves)
- self.playlist_generics={}
-
+ self.playlist_generics={}
+
#make sure we execute _plug_init() for every command line plugin we import
for plugin_name in self.config['plugins']:
try:
#load default list, if possible
self.do_loadlist(self.config['defaultlist'])
-
+
#HELPER FUNCTIONS
#Everything sending an event should be here
def _measure_N_points(self, N, whatset=1):
except Empty:
pass
return points
-
+
def _get_displayed_plot(self,dest=0):
'''
returns the currently displayed plot.
if displayed_plot:
break
return displayed_plot
-
+
def _send_plot(self,plots):
'''
sends a plot to the GUI
'''
wx.PostEvent(self.frame, self.list_of_events['plot_graph'](plots=plots))
return
-
+
def _find_plotmanip(self, name):
'''
returns a plot manipulator function from its name
'''
return self.plotmanip[self.config['plotmanips'].index(name)]
-
+
def _clickize(self, xvector, yvector, index):
'''
- returns a ClickedPoint() object from an index and vectors of x, y coordinates
+ returns a ClickedPoint() object from an index and vectors of x, y coordinates
'''
point=ClickedPoint()
point.index=index
point.absolute_coords=xvector[index],yvector[index]
point.find_graph_coords(xvector,yvector)
return point
-
+
#HERE COMMANDS BEGIN
-
+
def help_set(self):
print '''
SET
value=None
else:
value=args[1]
-
+
self.config[key]=value
self.do_plot(0)
-
+
#PLAYLIST MANAGEMENT AND NAVIGATION
#------------------------------------
-
+
def help_loadlist(self):
print '''
LOADLIST
#checking for args: if nothing is given as input, we warn and exit.
while len(args)==0:
args=linp.safeinput('File to load?')
-
+
arglist=args.split()
play_to_load=arglist[0]
-
+
#We assume a Hooke playlist has the extension .hkp
if play_to_load[-4:] != '.hkp':
play_to_load+='.hkp'
-
- try:
+
+ try:
playxml=PlaylistXML()
self.current_list, self.playlist_generics=playxml.load(play_to_load)
self.current_playxml=playxml
except IOError:
print 'File not found.'
return
-
+
print 'Loaded %s curves' %len(self.current_list)
-
+
if 'pointer' in self.playlist_generics.keys():
self.pointer=int(self.playlist_generics['pointer'])
else:
#if no pointer is found, set the current curve as the first curve of the loaded playlist
self.pointer=0
print 'Starting at curve ',self.pointer
-
+
self.current=self.current_list[self.pointer]
-
+
#resets saved/notes saved state
self.playlist_saved=0
self.playlist_name=''
- self.notes_saved=0
-
+ self.notes_saved=0
+
self.do_plot(0)
-
-
+
+
def help_genlist(self):
print '''
GENLIST
are all equivalent syntax.
------------
Syntax: genlist [input files]
-
+
'''
def do_genlist(self,args):
#args list is: input path, output name
if len(args)==0:
args=linp.safeinput('Input files?')
-
- arglist=args.split()
+
+ arglist=args.split()
list_path=arglist[0]
-
+
#if it's a directory, is like /directory/*.*
#FIXME: probably a bit kludgy.
- if os.path.isdir(list_path):
+ if os.path.isdir(list_path):
if platform.system == 'Windows':
SLASH="\\"
else:
SLASH="/"
if list_path[-1] == SLASH:
list_path=list_path+'*.*'
- else:
+ else:
list_path=list_path+SLASH+'*.*'
-
- #expanding correctly the input list with the glob module :)
+
+ #expanding correctly the input list with the glob module :)
list_files=glob.glob(list_path)
list_files.sort()
for item in list_files:
try:
if os.path.isfile(item):
- self.current_list.append(lhc.HookeCurve(os.path.abspath(item)))
+ self.current_list.append(lhc.HookeCurve(os.path.abspath(item)))
except:
pass
-
- self.pointer=0
+
+ self.pointer=0
if len(self.current_list)>0:
self.current=self.current_list[self.pointer]
else:
print 'Empty list!'
return
-
+
#resets saved/notes saved state
self.playlist_saved=0
self.playlist_name=''
- self.notes_saved=0
-
+ self.notes_saved=0
+
self.do_plot(0)
-
-
+
+
def do_savelist(self,args):
'''
SAVELIST
'''
while len(args)==0:
args=linp.safeinput('Output file?',['savedlist.txt'])
-
+
output_filename=args
-
+
self.playlist_generics['pointer']=self.pointer
-
+
#autocomplete filename if not specified
if output_filename[-4:] != '.hkp':
output_filename+='.hkp'
-
+
playxml=PlaylistXML()
playxml.export(self.current_list, self.playlist_generics)
- playxml.save(output_filename)
-
+ playxml.save(output_filename)
+
#remembers we have saved playlist
self.playlist_saved=1
-
+
def help_addtolist(self):
print '''
ADDTOLIST
print 'You must give the input filename you want to add'
self.help_addtolist()
return
-
+
filenames=glob.glob(args)
-
+
for filename in filenames:
self.current_list.append(lhc.HookeCurve(os.path.abspath(filename)))
#we need to save playlist
self.playlist_saved=0
-
+
def help_printlist(self):
print '''
PRINTLIST
def do_printlist(self,args):
for item in self.current_list:
print item.path
-
-
+
+
def help_jump(self):
print '''
JUMP
Syntax: jump {$curve}
If the curve is not in the current playlist, it politely asks if we want to add it.
- '''
+ '''
def do_jump(self,filename):
'''
jumps to the curve with the given filename.
if the filename is not in the playlist, it asks if we must add it or not.
'''
-
+
if filename=='':
filename=linp.safeinput('Jump to?')
-
+
filepath=os.path.abspath(filename)
print filepath
-
+
c=0
item_not_found=1
while item_not_found:
try:
-
+
if self.current_list[c].path == filepath:
self.pointer=c
self.current=self.current_list[self.pointer]
item_not_found=0
self.do_plot(0)
else:
- c+=1
+ c+=1
except IndexError:
#We've found the end of the list.
answer=linp.safeinput('Curve not found in playlist. Add it to list?',['y'])
self.current=self.current_list[-1]
self.pointer=(len(current_list)-1)
self.do_plot(0)
-
+
item_not_found=0
-
-
+
+
def do_index(self,args):
'''
INDEX
-----
Syntax: index
'''
- print self.pointer+1, 'of', len(self.current_list)
-
-
+ print self.pointer+1, 'of', len(self.current_list)
+
+
def help_next(self):
print '''
NEXT
print 'No curve file loaded, currently!'
print 'This should not happen, report to http://code.google.com/p/hooke'
return
-
+
if self.pointer == (len(self.current_list)-1):
self.pointer=0
print 'Playlist finished; back to first curve.'
else:
self.pointer+=1
-
+
self.current=self.current_list[self.pointer]
self.do_plot(0)
-
-
+
+
def help_n(self):
self.help_next()
def do_n(self,args):
self.do_next(args)
-
+
def help_previous(self,args):
print '''
PREVIOUS
return
if self.pointer == 0:
self.pointer=(len(self.current_list)-1)
- print 'Start of playlist; jump to last curve.'
+ print 'Start of playlist; jump to last curve.'
else:
self.pointer-=1
-
+
self.current=self.current_list[self.pointer]
self.do_plot(args)
-
-
+
+
def help_p(self):
self.help_previous()
def do_p(self,args):
self.do_previous(args)
-
+
#PLOT INTERACTION COMMANDS
-#-------------------------------
+#-------------------------------
def help_plot(self):
print '''
PLOT
Syntax: plot
'''
def do_plot(self,args):
-
+
self.current.identify(self.drivers)
self.plots=self.current.curve.default_plots()
try:
print 'Unexpected error occurred in do_plot().'
print e
return
-
+
#apply the plotmanip functions eventually present
nplots=len(self.plots)
c=0
while c<nplots:
for function in self.plotmanip: #FIXME: something strange happens about self.plotmanip[0]
self.plots[c]=function(self.plots[c], self.current)
-
+
self.plots[c].xaxes=self.config['xaxes'] #FIXME: in the future, xaxes and yaxes should be set per-plot
self.plots[c].yaxes=self.config['yaxes']
-
+
c+=1
self._send_plot(self.plots)
-
+
def _delta(self, set=1):
'''
calculates the difference between two clicked points
unitx=self.plots[points[0].dest].units[0]
unity=self.plots[points[0].dest].units[1]
return dx,unitx,dy,unity
-
+
def do_delta(self,args):
'''
DELTA
-
+
Measures the delta X and delta Y between two points.
----
Syntax: delta
dx,unitx,dy,unity=self._delta()
print str(dx)+' '+unitx
print str(dy)+' '+unity
-
+
def _point(self, set=1):
'''calculates the coordinates of a single clicked point'''
print 'Click one point'
point=self._measure_N_points(N=1, whatset=set)
-
+
x=point[0].graph_coords[0]
y=point[0].graph_coords[1]
unitx=self.plots[point[0].dest].units[0]
unity=self.plots[point[0].dest].units[1]
return x,unitx,y,unity
-
+
def do_point(self,args):
'''
POINT
-
+
Returns the coordinates of a point on the graph.
----
Syntax: point
print str(x)+' '+unitx
print str(y)+' '+unity
to_dump='point '+self.current.path+' '+str(x)+' '+unitx+', '+str(y)+' '+unity
- self.outlet.push(to_dump)
-
-
+ self.outlet.push(to_dump)
+
+
def do_close(self,args=None):
'''
CLOSE
to_close=1
else:
to_close=1
-
+
close_plot=self.list_of_events['close_plot']
wx.PostEvent(self.frame, close_plot(to_close=to_close))
-
+
def do_show(self,args=None):
'''
SHOW
Shows both plots.
- '''
+ '''
show_plots=self.list_of_events['show_plots']
wx.PostEvent(self.frame, show_plots())
-
-
-
+
+
+
#PLOT EXPORT AND MANIPULATION COMMANDS
def help_export(self):
print '''
'''
def do_export(self,args):
#FIXME: the bottom plot doesn't have the title
-
+
dest=0
-
+
if len(args)==0:
#FIXME: We have to go into the libinput stuff and fix it, for now here's a dummy replacement...
#name=linp.safeinput('Filename?',[self.current.path+'.png'])
args=args.split()
name=args[0]
if len(args) > 1:
- dest=int(args[1])
-
+ dest=int(args[1])
+
export_image=self.list_of_events['export_image']
wx.PostEvent(self.frame, export_image(name=name, dest=dest))
-
-
+
+
def help_txt(self):
print '''
TXT
Syntax: txt [filename] {plot to export}
'''
def do_txt(self,args):
-
+
def transposed2(lists, defval=0):
'''
transposes a list of lists, i.e. from [[a,b,c],[x,y,z]] to [[a,x],[b,y],[c,z]] without losing
'''
if not lists: return []
return map(lambda *row: [elem or defval for elem in row], *lists)
-
+
whichplot=0
args=args.split()
if len(args)==0:
whichplot=int(args[1])
except:
pass
-
- columns=[]
+
+ columns=[]
for dataset in self.plots[whichplot].vectors:
- for i in range(0,len(dataset)):
+ for i in range(0,len(dataset)):
columns.append([])
for value in dataset[i]:
- columns[-1].append(str(value))
-
+ columns[-1].append(str(value))
+
rows=transposed2(columns, 'nan')
rows=[' , '.join(item) for item in rows]
text='\n'.join(rows)
-
+
txtfile=open(filename,'w+')
#Save units of measure in header
txtfile.write('X:'+self.plots[whichplot].units[0]+'\n')
txtfile.write('Y:'+self.plots[whichplot].units[1]+'\n')
txtfile.write(text)
txtfile.close()
-
-
+
+
#LOGGING, REPORTING, NOTETAKING
-
+
def do_note_old(self,args):
'''
NOTE_OLD
**deprecated**: Use note instead. Will be removed in 0.9
-
+
Writes or displays a note about the current curve.
If [anything] is empty, it displays the note, otherwise it adds a note.
The note is then saved in the playlist if you issue a savelist command
---------------
- Syntax: note_old [anything]
+ Syntax: note_old [anything]
'''
if args=='':
args=args.decode('ascii','ignore')
if len(args)==0:
args='?'
-
+
self.current_list[self.pointer].notes=args
self.notes_saved=0
-
-
+
+
def do_note(self,args):
'''
NOTE
-
+
Writes or displays a note about the current curve.
If [anything] is empty, it displays the note, otherwise it adds a note.
The note is then saved in the playlist if you issue a savelist command.
---------------
- Syntax: note_old [anything]
+ Syntax: note_old [anything]
'''
if args=='':
f=open(self.notes_filename,'w')
f.write(title_line)
f.close()
-
- #bypass UnicodeDecodeError troubles
+
+ #bypass UnicodeDecodeError troubles
try:
args=args.decode('ascii')
except:
if len(args)==0:
args='?'
self.current_list[self.pointer].notes=args
-
+
f=open(self.notes_filename,'a+')
note_string=(self.current.path+' | '+self.current.notes+'\n')
f.write(note_string)
f.close()
-
+
def help_notelog(self):
print '''
NOTELOG
Writes a log of the notes taken during the session for the current
playlist
---------------
+--------------
Syntax notelog [filename]
-'''
+'''
def do_notelog(self,args):
-
+
if len(args)==0:
args=linp.safeinput('Notelog filename?',['notelog.txt'])
-
+
note_lines='Notes taken at '+time.asctime()+'\n'
for item in self.current_list:
if len(item.notes)>0:
#FIXME: file path should be truncated...
note_string=(item.path+' | '+item.notes+'\n')
note_lines+=note_string
-
+
try:
f=open(args,'a+')
f.write(note_lines)
except IOError, (ErrorNumber, ErrorMessage):
print 'Error: notes cannot be saved. Catched exception:'
print ErrorMessage
-
+
self.notes_saved=1
def help_copylog(self):
Syntax copylog [directory]
'''
def do_copylog(self,args):
-
+
if len(args)==0:
args=linp.safeinput('Destination directory?') #TODO default
-
+
mydir=os.path.abspath(args)
if not os.path.isdir(mydir):
print 'Destination is not a directory.'
return
-
+
for item in self.current_list:
if len(item.notes)>0:
try:
self.outlet.delete(args)
#OS INTERACTION COMMANDS
-#-----------------
+#-----------------
def help_dir(self):
print '''
DIR, LS
ls [path]
'''
def do_dir(self,args):
-
+
if len(args)==0:
args='*'
print glob.glob(args)
-
+
def help_ls(self):
self.help_dir(self)
def do_ls(self,args):
self.do_dir(args)
-
+
def help_pwd(self):
print '''
PWD
Syntax: pwd
'''
def do_pwd(self,args):
- print os.getcwd()
-
+ print os.getcwd()
+
def help_cd(self):
print '''
CD
os.chdir(mypath)
except OSError:
print 'I cannot access that directory.'
-
-
+
+
def help_system(self):
print '''
SYSTEM
'''
pass
def do_system(self,args):
- waste=os.system(args)
-
+ waste=os.system(args)
+
def do_debug(self,args):
'''
this is a dummy command where I put debugging things
'''
print self.config['plotmanips']
pass
-
+
def help_current(self):
print '''
CURRENT
'''
def do_current(self,args):
print self.current.path
-
+
def do_info(self,args):
'''
INFO
for set in plot.vectors:
lengths=[len(item) for item in set]
print 'Data set size: ',lengths
-
+
def do_version(self,args):
'''
VERSION
------
Prints the current version and codename, plus library version. Useful for debugging.
- '''
+ '''
print 'Hooke '+__version__+' ('+__codename__+')'
print 'Released on: '+__releasedate__
print '---'
print 'Platform: '+str(platform.uname())
print '---'
print 'Loaded plugins:',self.config['loaded_plugins']
-
+
def help_exit(self):
print '''
EXIT, QUIT
------
Syntax: exit
Syntax: quit
-'''
+'''
def do_exit(self,args):
we_exit='N'
-
+
if (not self.playlist_saved) or (not self.notes_saved):
we_exit=linp.safeinput('You did not save your playlist and/or notes. Exit?',['n'])
else:
we_exit=linp.safeinput('Exit?',['y'])
-
+
if we_exit[0].upper()=='Y':
wx.CallAfter(self.frame.Close)
sys.exit(0)
else:
return
-
+
def help_quit(self):
self.help_exit()
def do_quit(self,args):
self.do_exit(args)
-
-
-
if __name__ == '__main__':
mycli=HookeCli(0)
mycli.cmdloop()
#!/usr/bin/env python
-
'''
libhooke.py
This program is released under the GNU General Public License version 2.
'''
-
-
import libhookecurve as lhc
import scipy
import csv
HOOKE_VERSION=['0.8.3_devel', 'Seinei', '2008-04-16']
-WX_GOOD=['2.6','2.8']
-
+WX_GOOD=['2.6','2.8']
+
class PlaylistXML:
'''
This module allows for import/export of an XML playlist into/out of a list of HookeCurve objects
'''
-
+
def __init__(self):
-
+
self.playlist=None #the DOM object representing the playlist data structure
self.playpath=None #the path of the playlist XML file
self.plaything=None
self.hidden_attributes=['curve'] #This list contains hidden attributes that we don't want to go into the playlist.
-
+
def export(self, list_of_hooke_curves, generics):
'''
Creates an initial playlist from a list of files.
<element path="/my/file/path/"/ attribute="attribute">
<element path="...">
</playlist>
- '''
-
+ '''
+
#create the output playlist, a simple XML document
impl=xml.dom.minidom.getDOMImplementation()
#create the document DOM object and the root element
newdoc=impl.createDocument(None, "playlist",None)
top_element=newdoc.documentElement
-
+
#save generics variables
playlist_generics=newdoc.createElement("generics")
top_element.appendChild(playlist_generics)
for key in generics.keys():
newdoc.createAttribute(key)
playlist_generics.setAttribute(key,str(generics[key]))
-
+
#save curves and their attributes
for item in list_of_hooke_curves:
#playlist_element=newdoc.createElement("curve")
for key in item.__dict__:
if not (key in self.hidden_attributes):
newdoc.createAttribute(key)
- playlist_element.setAttribute(key,str(item.__dict__[key]))
-
+ playlist_element.setAttribute(key,str(item.__dict__[key]))
+
self.playlist=newdoc
-
+
def load(self,filename):
'''
loads a playlist file
'''
myplay=file(filename)
self.playpath=filename
-
+
#the following 3 lines are needed to strip newlines. otherwise, since newlines
#are XML elements too (why?), the parser would read them (and re-save them, multiplying
#newlines...)
the_file=myplay.read()
the_file_lines=the_file.split('\n')
the_file=''.join(the_file_lines)
-
- self.playlist=xml.dom.minidom.parseString(the_file)
-
+
+ self.playlist=xml.dom.minidom.parseString(the_file)
+
#inner parsing functions
def handlePlaylist(playlist):
list_of_files=playlist.getElementsByTagName("element")
generics=playlist.getElementsByTagName("generics")
return handleFiles(list_of_files), handleGenerics(generics)
-
+
def handleGenerics(generics):
generics_dict={}
if len(generics)==0:
return generics_dict
-
+
for attribute in generics[0].attributes.keys():
generics_dict[attribute]=generics[0].getAttribute(attribute)
return generics_dict
-
+
def handleFiles(list_of_files):
new_playlist=[]
for myfile in list_of_files:
for attribute in myfile.attributes.keys(): #extract attributes for the single curve
the_curve.__dict__[attribute]=myfile.getAttribute(attribute)
new_playlist.append(the_curve)
-
+
return new_playlist #this is the true thing returned at the end of this function...(FIXME: clarity)
-
+
return handlePlaylist(self.playlist)
-
+
def save(self,output_filename):
'''
except IOError:
print 'libhooke.py : Cannot save playlist. Wrong path or filename'
return
-
+
self.playlist.writexml(outfile,indent='\n')
outfile.close()
class HookeConfig:
'''
Handling of Hooke configuration file
-
+
Mostly based on the simple-yet-useful examples of the Python Library Reference
about xml.dom.minidom
-
+
FIXME: starting to look a mess, should require refactoring
'''
-
+
def __init__(self):
self.config={}
self.config['plugins']=[]
self.config['drivers']=[]
self.config['plotmanips']=[]
-
+
def load_config(self, filename):
myconfig=file(filename)
-
+
#the following 3 lines are needed to strip newlines. otherwise, since newlines
#are XML elements too, the parser would read them (and re-save them, multiplying
#newlines...)
the_file=myconfig.read()
the_file_lines=the_file.split('\n')
the_file=''.join(the_file_lines)
-
- self.config_tree=xml.dom.minidom.parseString(the_file)
-
+
+ self.config_tree=xml.dom.minidom.parseString(the_file)
+
def getText(nodelist):
#take the text from a nodelist
#from Python Library Reference 13.7.2
if node.nodeType == node.TEXT_NODE:
rc += node.data
return rc
-
+
def handleConfig(config):
display_elements=config.getElementsByTagName("display")
plugins_elements=config.getElementsByTagName("plugins")
handleWorkdir(workdir_elements)
handleDefaultlist(defaultlist_elements)
handlePlotmanip(plotmanip_elements)
-
+
def handleDisplay(display_elements):
for element in display_elements:
for attribute in element.attributes.keys():
self.config[attribute]=element.getAttribute(attribute)
-
+
def handlePlugins(plugins):
for plugin in plugins[0].childNodes:
try:
self.config['drivers'].append(str(driver.tagName))
except: #if we allow fancy formatting of xml, there is a text node, so tagName fails for it...
pass
-
+
def handlePlotmanip(plotmanips):
for plotmanip in plotmanips[0].childNodes:
try:
self.config['plotmanips'].append(str(plotmanip.tagName))
except: #if we allow fancy formatting of xml, there is a text node, so tagName fails for it...
pass
-
+
def handleWorkdir(workdir):
'''
default working directory
'''
wdir=getText(workdir[0].childNodes)
self.config['workdir']=wdir.strip()
-
+
def handleDefaultlist(defaultlist):
'''
default playlist
'''
dflist=getText(defaultlist[0].childNodes)
self.config['defaultlist']=dflist.strip()
-
+
handleConfig(self.config_tree)
#making items in the dictionary more machine-readable
for item in self.config.keys():
self.config[item]=None
else:
pass
-
+
return self.config
-
-
+
+
def save_config(self, config_filename):
print 'Not Implemented.'
- pass
+ pass
class ClickedPoint:
this class defines what a clicked point on the curve plot is
'''
def __init__(self):
-
+
self.is_marker=None #boolean ; decides if it is a marker
self.is_line_edge=None #boolean ; decides if it is the edge of a line (unused)
self.absolute_coords=(None,None) #(float,float) ; the absolute coordinates of the clicked point on the graph
self.graph_coords=(None,None) #(float,float) ; the coordinates of the plot that are nearest in X to the clicked point
self.index=None #integer ; the index of the clicked point with respect to the vector selected
self.dest=None #0 or 1 ; 0=top plot 1=bottom plot
-
-
+
+
def find_graph_coords_old(self, xvector, yvector):
'''
Given a clicked point on the plot, finds the nearest point in the dataset (in X) that
corresponds to the clicked point.
OLD & DEPRECATED - to be removed
'''
-
+
#FIXME: a general algorithm using min() is needed!
print '---DEPRECATED FIND_GRAPH_COORDS_OLD---'
best_index=0
best_dist=10**9 #should be more than enough given the scale
-
+
for index in scipy.arange(1,len(xvector),1):
dist=((self.absolute_coords[0]-xvector[index])**2)+(100*((self.absolute_coords[1]-yvector[index])))**2
#TODO, generalize? y coordinate is multiplied by 100 due to scale differences in the plot
if dist<best_dist:
best_index=index
best_dist=dist
-
+
self.index=best_index
self.graph_coords=(xvector[best_index],yvector[best_index])
return
-
+
def find_graph_coords(self,xvector,yvector):
'''
Given a clicked point on the plot, finds the nearest point in the dataset that
dists=[]
for index in scipy.arange(1,len(xvector),1):
dists.append(((self.absolute_coords[0]-xvector[index])**2)+((self.absolute_coords[1]-yvector[index])**2))
-
+
self.index=dists.index(min(dists))
self.graph_coords=(xvector[self.index],yvector[self.index])
#-----------------------------------------
-#CSV-HELPING FUNCTIONS
-
+#CSV-HELPING FUNCTIONS
+
def transposed2(lists, defval=0):
'''
transposes a list of lists, i.e. from [[a,b,c],[x,y,z]] to [[a,x],[b,y],[c,z]] without losing
'''
if not lists: return []
return map(lambda *row: [elem or defval for elem in row], *lists)
-
+
def csv_write_dictionary(f, data, sorting='COLUMNS'):
'''
Writes a CSV file from a dictionary, with keys as first column or row
Keys are in "random" order.
-
+
Keys should be strings
Values should be lists or other iterables
'''
writer.writerow(keys)
for item in t_values:
writer.writerow(item)
-
+
if sorting=='ROWS':
print 'Not implemented!' #FIXME: implement it.
-#-----------------------------------------
-
+#-----------------------------------------
+
def debug():
'''
debug stuff from latest rewrite of hooke_playlist.py
print confo.load_config('hooke.conf')
if __name__ == '__main__':
- debug()
\ No newline at end of file
+ debug()
-#!/usr/bin/env python
-
-
class HookeCurve(object):
-
+
def __init__(self,path):
self.path=path
self.curve=Driver()
self.notes=''
-
+
def identify(self, drivers):
'''
identifies a curve and returns the corresponding object
self.curve=tempcurve
del tempcurve
return True
-
+
print 'Not a recognizable curve format.'
return False
-
-
+
+
class Driver:
'''
Base class for file format drivers.
-
+
To be overridden
'''
def __init__(self):
self.experiment=''
self.filetype=''
-
+
def is_me(self):
'''
This method must read the file and return True if the filetype can be managed by the driver, False if not.
'''
return False
-
+
def close_all(self):
'''
This method must close all the open files of the driver, explicitly.
'''
return None
-
+
def default_plots(self):
dummy_default=PlotObject()
dummy_default.vectors.append([[[0]],[[0]]])
return [dummy_default]
-
+
class PlotObject:
-
+
def __init__(self):
-
+
'''
the plot destination
0=top
1=bottom
'''
- self.destination=0
-
+ self.destination=0
+
'''
self.vectors is a multidimensional array:
self.vectors[0]=plot1
self.vectors[1]=plot2
self.vectors[2]=plot3
etc.
-
+
2 curves in a x,y plot are:
[[[x1],[y1]],[[x2],[y2]]]
for example:
self.units is simpler. for each plot with N axes (x,y,z...) only N labels
can be made, regardless of the number of superimposed plots
so units for the double plot above is: [unitx, unity]
-
+
units are strings
'''
self.units=['','']
-
+
'''
xaxes and yaxes directions. 0,0 means the common +X=right, +Y=top directions
'''
self.xaxes=0
self.yaxes=0
-
- self.title='' #title
-
+
+ self.title='' #title
+
'''
styles: defines what is the style of the current plots. If undefined or None, it is line plot.
If an element of the list is 'scatter', the corresponding dataset
is drawn with scattered points and not a continuous line.
'''
self.styles=[]
-
+
'''
colors: define what is the colour of the current plots
'''
self.colors=[]
-
+
def add_set(self,x,y):
'''
Adds an x,y data set to the vectors.
self.vectors[-1].append(x)
self.vectors[-1].append(y)
return
-
+
def remove_set(self,whichset):
'''
Removes a set
'''
waste=self.vectors.pop(whichset)
return
-
+
def normalize_vectors(self):
'''
Trims the vector lengths as to be equal in a plot.
'''
-
+
for index in range(0,len(self.vectors)):
vectors_to_plot=self.vectors[index]
lengths=[len(vector) for vector in vectors_to_plot]
-#!/usr/bin/env python
-
'''
Input check routines.
from types import *
-
def safeinput (message, valid=[]):
'''
friendlier frontend for alphainput and numinput
valid should be a list of 0...n values
'''
- #if possible values are not listed we just ask for any non-null input
+ #if possible values are not listed we just ask for any non-null input
if len(valid)==0:
return alphainput(message, '',1,[])
-
-
+
+
if len(valid)>0:
#if valid values are string we use alphainput, if it is only one we take as default
if type(valid[0]) is StringType:
return alphainput(message, valid[0], 0,[])
else:
return alphainput(message,'', 1,valid)
-
+
#if valid values are numbers we use numinput
if type(valid[0]) is IntType:
if len(valid)==1:
return numinput(message,valid[0],1,[])
else:
return numinput(message,'',1,valid)
-
-
+
def alphainput (message, default, repeat, valid):
'''
default: return value if user input was not correct (and repeat=0)
repeat: keeps asking user till it gets a valid input
valid: list of allowed answers, empty list for "anything"
- '''
+ '''
if default and not repeat:
print 'Press [enter] for default: ('+str(default)+')'
reply=raw_input(message)
if len(valid)>0:
- if reply in valid:
+ if reply in valid:
return reply
else:
if repeat==1:
reply=raw_input(message)
return reply
-
def checkalphainput (test, default, valid):
#useful when input was taken form command args
if len(valid)>0:
- if test in valid:
+ if test in valid:
return test
else:
return default
default: return value if user input was not correct (and repeat=0)
repeat: keeps asking user till it gets a valid input
limits: pair of values, input is checked to be between them, empty list for "any number"
- '''
+ '''
if default and not repeat:
print 'Press [enter] for default: '+str(default)
-
+
reply=raw_input(message)
-
+
try:
intreply=int(reply)
except:
intreply=None
-
+
if len(limits)==2:
high=int(limits.pop())
low=int(limits.pop())
intreply=None
return intreply
+
def checknuminput(test,default,limits):
#useful when input was taken from command args
if len(limits)==2:
return int(test)
else:
return default
-
-#!/usr/bin/env python
-
'''
Basic outlet object
This program is released under the GNU General Public License version 2.
'''
-
import re
class Outlet(object):
-
+
def __init__(self):
self.buffer=[]
#relations is still unused
self.relations=[]
-
+
def push(self, args):
#queue new entry
self.buffer.append(args)
def empty(self):
self.buffer=[]
-
+
def read_last(self):
return self.buffer[len(self.buffer)-1]
-
+
def read_first(self):
return self.buffer[0]
def read_type(self,dtype):
- #returns entries matching a given type (force, distance, point...)
+ #returns entries matching a given type (force, distance, point...)
aux=[]
index=0
if dtype=='all':
if re.match(dtype+'*',i):
aux.append(i)
return aux
-
-
-#!/usr/bin/env python
-
'''
a library of helping functions for spotting force spectroscopy peaks
dim=len(data)
window=len(vect)
temparr=[0.0]*dim
-
+
end=dim-window
for j in range(end):
for k in range(window):
temparr[j]+=data[j+k]*vect[k]
- return temparr
-
+ return temparr
+
def absdev(arr):
'''
Calculates the absolute deviation of a vector
'''
med=0.0
absD=0.0
-
+
med=np.mean(arr)
for j in arr:
absD+=abs(j-med)
return absD/len(arr)
-
+
def noise_absdev(data,positive=False,maxcut=0.2,stable=0.005):
'''
Returns the standard deviation of the noise.
The logic is: we cut the most negative (or positive) data points until the absolute deviation
becomes stable (it doesn't vary more than 0.005) or we have cut more than maxcut*len(data) points.
Then calculate the absolute deviation.
-
+
If positive=True we cut the most positive data points, False=we cut the negative ones.
'''
out=[item for item in data] #we copy just to be sure...
out.sort()
if positive:
out.reverse()
-
+
temp_absdev=absdev(out)
-
+
for index in range(len(out)):
cutindex=(index+1)*5
cut_absdev=absdev(out[cutindex:]) #we jump five points after five points...
temp_absdev=cut_absdev
else:
break
-
+
return cut_absdev
-
+
def abovenoise(convoluted,noise_level,cut_index=0,abs_devs=4):
'''
Generates a vector which is 0 where the vector is less than abs_devs*noise_level ; 1 if not (spike).
'''
#calculate absolute noise deviation
#noise_level=noise_absdev(convoluted[cut_index:])
-
+
above=[]
-
+
for index in range(len(convoluted)):
if index<cut_index:
above.append(0)
else:
above.append(0)
return above
-
+
def find_peaks(above, seedouble=10):
'''
Finds individual peak location.
abovenoise() finds all points above a given threshold in the convolution. This point is often not unique
but there is a "cluster" of points above a threshold.
- Here we obtain only the coordinates of the largest *convolution* spike for each cluster.
-
+ Here we obtain only the coordinates of the largest *convolution* spike for each cluster.
+
above=vector obtained from abovenoise()
seedouble=value at which we want to "delete" double peaks. That is, if two peaks have a distance
< than $seedouble points , only the first is kept.
nonzero=[]
peaks_location=[]
peaks_size=[]
-
+
for index in range(len(above)):
if above[index] != 0:
nonzero.append(index)
nonzero=[]
else:
pass
-
+
#recursively eliminate double peaks
#(maybe not the smartest of cycles, but i'm asleep...)
temp_location=None
if peaks_location[index+1]-peaks_location[index] < seedouble:
temp_location=peaks_location[:index]+peaks_location[index+1:]
if temp_location != []:
- peaks_location=temp_location
-
- return peaks_location,peaks_size
\ No newline at end of file
+ peaks_location=temp_location
+
+ return peaks_location,peaks_size
-#!/usr/bin/env python
-
'''
Basic Viewer and ascii saver example
This program is released under the GNU General Public License version 2.
'''
-
import liboutlet as lout
import libinput as linput
+
class Viewer(object):
source=[]
data=[]
dtype='all'
action=[] #alias to call the actual viewer function, makes it general
-
+
def setdtype(self, dt):
#determines wich type of data will be retrieved from outlet
self.data=self.source.read_type(self.dtype)
-
class Ascii(Viewer):
#example viewer, it just retrieves data and writes it to a text file
#probably this should be in viewer.py?
destfile=open(destination,'w+')
destfile.write('\n'.join(self.data))
destfile.close()
-
-#!/usr/bin/env python
# -*- coding: utf-8 -*-
from libhooke import WX_GOOD, ClickedPoint
class autopeakCommands:
-
+
def do_autopeak(self,args):
'''
AUTOPEAK
- measures peak maximum forces with a baseline
- measures slope in proximity of peak maximum
Requires flatten plotmanipulator , fit.py plugin , flatfilts.py plugin with convfilt
-
+
Syntax:
autopeak [rebase] [pl=value] [t=value] [noauto] [reclick]
-
+
rebase : Re-asks baseline interval
-
- pl=[value] : Use a fixed persistent length for the fit. If pl is not given,
- the fit will be a 2-variable
+
+ pl=[value] : Use a fixed persistent length 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 meters.
+ The value must be in meters.
Scientific notation like 0.35e-9 is fine.
-
+
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.
-
- noauto : allows for clicking the contact point by
+
+ noauto : allows for clicking the contact point by
hand (otherwise it is automatically estimated) the first time.
If subsequent measurements are made, the same contact point
clicked the first time is used
-
+
reclick : redefines by hand the contact point, if noauto has been used before
but the user is unsatisfied of the previously choosen contact point.
-
+
usepoints : fit interval by number of points instead than by nanometers
-
+
noflatten : does not use the "flatten" plot manipulator
-
+
When you first issue the command, it will ask for the filename. If you are giving the filename
of an existing file, autopeak will resume it and append measurements to it. If you are giving
a new filename, it will create the file and append to it until you close Hooke.
-
-
+
+
Useful variables (to set with SET command):
---
fit_function = type of function to use for elasticity. If "wlc" worm-like chain is used, if "fjc" freely jointed
chain is used
-
+
temperature= temperature of the system for wlc/fjc fit (in K)
-
+
auto_slope_span = number of points on which measure the slope, for slope
-
+
auto_fit_nm = number of nm to fit before the peak maximum, for WLC/FJC (if usepoints false)
auto_fit_points = number of points to fit before the peak maximum, for WLC/FJC (if usepoints true)
-
+
baseline_clicks = -1: no baseline, f=0 at the contact point (whether hand-picked or automatically found)
0: automatic baseline
1: decide baseline with a single click and length defined in auto_left_baseline
2: let user click points of baseline
auto_left_baseline = length in nm to use as baseline from the right point (if baseline_clicks=0 , 1)
auto_right_baseline = distance in nm of peak-most baseline point from last peak (if baseline_clicks = 0)
-
+
auto_min_p ; auto_max_p = Minimum and maximum persistence length (if using WLC) or Kuhn length (if using FJC)
outside of which the peak is automatically discarded (in nm)
'''
-
+
#MACROS.
#FIXME: to move outside function
def fit_interval_nm(start_index,plot,nm,backwards):
'''
whatset=1 #FIXME: should be decidable
x_vect=plot.vectors[1][0]
-
+
c=0
i=start_index
start=x_vect[start_index]
while abs(x_vect[i]-x_vect[start_index])*(10**9) < nm:
if i==0 or i==maxlen-1: #we reached boundaries of vector!
return c
-
+
if backwards:
i-=1
else:
i+=1
c+=1
return c
-
+
def pickup_contact_point():
'''macro to pick up the contact point by clicking'''
contact_point=self._measure_N_points(N=1, whatset=1)[0]
self.wlccontact_index=contact_point.index
self.wlccurrent=self.current.path
return contact_point, contact_point_index
-
+
def find_current_peaks(noflatten):
#Find peaks.
defplot=self.current.curve.default_plots()[0]
defplot=flatten(defplot, self.current, customvalue=1) #Flatten curve before feeding it to has_peaks
peak_location,peak_size=self.has_peaks(defplot, self.convfilt_config['mindeviation'])
return peak_location, peak_size
-
+
#default fit etc. variables
pl_value=None
T=self.config['temperature']
-
+
slope_span=int(self.config['auto_slope_span'])
delta_force=10
rebase=False #if true=we select rebase
noflatten=False #if true=we avoid flattening
-
+
#initialize output data vectors
c_lengths=[]
p_lengths=[]
sigma_p_lengths=[]
forces=[]
slopes=[]
-
+
#pick up plot
displayed_plot=self._get_displayed_plot(0)
-
+
#COMMAND LINE PARSING
#--Using points instead of nm interval
if 'usepoints' in args.split():
usepoints=False
#--Recalculate baseline
if 'rebase' in args or (self.basecurrent != self.current.path):
- rebase=True
-
+ rebase=True
+
if 'noflatten' in args:
noflatten=True
-
+
#--Custom persistent length / custom temperature
for arg in args.split():
#look for a persistent length argument.
#look for a T argument. FIXME: spaces are not allowed between 'pl' and value
if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
t_expression=arg.split('=')
- T=float(t_expression[1])
+ T=float(t_expression[1])
#--Contact point arguments
if 'reclick' in args.split():
print 'Click contact point'
cindex=self.find_contact_point()
contact_point=self._clickize(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1], cindex)
#--END COMMAND LINE PARSING--
-
-
+
+
peak_location, peak_size = find_current_peaks(noflatten)
-
+
if len(peak_location) == 0:
print 'No peaks to fit.'
return
-
+
fitplot=copy.deepcopy(displayed_plot)
-
+
#Pick up force baseline
if rebase:
clicks=self.config['baseline_clicks']
self.basepoints.append(self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],base_index_1))
else:
self.basepoints=self._measure_N_points(N=2, whatset=whatset)
-
+
self.basecurrent=self.current.path
-
+
boundaries=[self.basepoints[0].index, self.basepoints[1].index]
boundaries.sort()
to_average=displayed_plot.vectors[1][1][boundaries[0]:boundaries[1]] #y points to average
avg=np.mean(to_average)
-
+
clicks=self.config['baseline_clicks']
if clicks==-1:
try:
avg=displayed_plot.vectors[1][1][contact_point_index]
except:
avg=displayed_plot.vectors[1][1][cindex]
-
+
for peak in peak_location:
#WLC FITTING
#define fit interval
fit_points=fit_interval_nm(peak, displayed_plot, self.config['auto_fit_nm'], True)
peak_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak)
other_fit_point=self._clickize(displayed_plot.vectors[1][0],displayed_plot.vectors[1][1],peak-fit_points)
-
+
#points for the fit
points=[contact_point, peak_point, other_fit_point]
-
+
if abs(peak_point.index-other_fit_point.index) < 2:
continue
-
+
if self.config['fit_function']=='wlc':
-
+
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)
elif self.config['fit_function']=='fjc':
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)
print 'Unknown fit function'
print 'Please set fit_function as wlc or fjc'
return
-
-
+
+
#Measure forces
delta_to_measure=displayed_plot.vectors[1][1][peak-delta_force:peak+delta_force]
y=min(delta_to_measure)
- #save force values (pN)
+ #save force values (pN)
#Measure slopes
slope=self.linefit_between(peak-slope_span,peak)[0]
-
-
+
+
#check fitted data and, if right, add peak to the measurement
#FIXME: code duplication
if len(params)==1: #if we did choose 1-value fit
sigma_p_lengths.append(0)
sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
forces.append(abs(y-avg)*(1.0e+12))
- slopes.append(slope)
+ slopes.append(slope)
#Add WLC fit lines to plot
fitplot.add_set(xfit,yfit)
if len(fitplot.styles)==0:
p_leng=params[1]*(1.0e+9)
#check if persistent length makes sense. otherwise, discard peak.
if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
- p_lengths.append(p_leng)
+ p_lengths.append(p_leng)
c_lengths.append(params[0]*(1.0e+9))
sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
forces.append(abs(y-avg)*(1.0e+12))
- slopes.append(slope)
-
+ slopes.append(slope)
+
#Add WLC fit lines to plot
fitplot.add_set(xfit,yfit)
if len(fitplot.styles)==0:
fitplot.colors.append(None)
else:
pass
-
-
+
+
#add basepoints to fitplot
- 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]])
+ 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]])
fitplot.styles.append('scatter')
fitplot.colors.append(None)
-
+
#Show wlc fits and peak locations
self._send_plot([fitplot])
#self.do_peaks('')
-
+
print 'Using fit function: ',self.config['fit_function']
print 'Measurements for all peaks detected:'
print 'contour (nm)', c_lengths
print 'sigma p (nm)',sigma_p_lengths
print 'forces (pN)',forces
print 'slopes (N/m)',slopes
-
+
controller=False
while controller==False:
#Ask the user what peaks to ignore from analysis.
print 'Bad input.'
controller=False
- #Clean data vectors from ignored peaks
+ #Clean data vectors from ignored peaks
#FIXME:code duplication
c_lengths=[item for item in c_lengths if item != None]
p_lengths=[item for item in p_lengths if item != None]
forces=[item for item in forces if item != None]
- slopes=[item for item in slopes if item != None]
- sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
- sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
-
+ slopes=[item for item in slopes if item != None]
+ sigma_c_lengths=[item for item in sigma_c_lengths if item != None]
+ sigma_p_lengths=[item for item in sigma_p_lengths if item != None]
+
print 'Measurements for chosen peaks:'
print 'contour (nm)',c_lengths
print 'sigma contour (nm)',sigma_c_lengths
print 'sigma p (nm)',sigma_p_lengths
print 'forces (pN)',forces
print 'slopes (N/m)',slopes
-
+
#Save file info
if self.autofile=='':
self.autofile=raw_input('Autopeak filename? (return to ignore) ')
if self.autofile=='':
print 'Not saved.'
return
-
+
if not os.path.exists(self.autofile):
f=open(self.autofile,'w+')
f.write('Analysis started '+time.asctime()+'\n')
f.write('----------------------------------------\n')
f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) ; Sigma contour (nm) ; Sigma persistence (nm)\n')
f.close()
-
+
print 'Saving...'
f=open(self.autofile,'a+')
-
+
f.write(self.current.path+'\n')
for i in range(len(c_lengths)):
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')
-
+
f.close()
self.do_note('autopeak')
-
\ No newline at end of file
for i in range(len(yarr)):
f.write(str(xarr[i])+";"+str(yarr[i])+"\n")
f.close()
-
-
-
-#!/usr/bin/env python
-
'''
FIT
class fitCommands:
-
+
def _plug_init(self):
self.wlccurrent=None
self.wlccontact_point=None
self.wlccontact_index=None
-
+
def wlc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
'''
Worm-like chain model fitting.
The function is the simple polynomial worm-like chain as proposed by C.Bustamante, J.F.Marko, E.D.Siggia
and S.Smith (Science. 1994 Sep 9;265(5178):1599-600.)
'''
-
+
'''
clicked_points[0] = contact point (calculated or hand-clicked)
clicked_points[1] and [2] are edges of chunk
'''
-
+
#STEP 1: Prepare the vectors to apply the fit.
-
-
+
+
if pl_value is not None:
pl_value=pl_value/(10**9)
-
+
#indexes of the selected chunk
first_index=min(clicked_points[1].index, clicked_points[2].index)
last_index=max(clicked_points[1].index, clicked_points[2].index)
-
+
#getting the chunk and reverting it
xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
xchunk.reverse()
- ychunk.reverse()
+ ychunk.reverse()
#put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
-
+
#make them arrays
xchunk_corr_up=scipy.array(xchunk_corr_up)
ychunk_corr_up=scipy.array(ychunk_corr_up)
-
-
+
+
#STEP 2: actually do the fit
-
+
#Find furthest point of chunk and add it a bit; the fit must converge
#from an excess!
xchunk_high=max(xchunk_corr_up)
xchunk_high+=(xchunk_high/10)
-
+
#Here are the linearized start parameters for the WLC.
#[lambd=1/Lo , pii=1/P]
-
+
p0=[(1/xchunk_high),(1/(3.5e-10))]
p0_plfix=[(1/xchunk_high)]
'''
ODR STUFF
fixme: remove these comments after testing
'''
-
-
+
+
def f_wlc(params,x,T=T):
'''
wlc function for ODR fitting
therm=Kb*T
y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
return y
-
+
def f_wlc_plfix(params,x,pl_value=pl_value,T=T):
'''
wlc function for ODR fitting
therm=Kb*T
y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
return y
-
+
#make the ODR fit
realdata=scipy.odr.RealData(xchunk_corr_up,ychunk_corr_up)
if pl_value:
else:
model=scipy.odr.Model(f_wlc)
o = scipy.odr.ODR(realdata, model, p0)
-
+
o.set_job(fit_type=2)
out=o.run()
fit_out=[(1/i) for i in out.beta]
-
+
#Calculate fit errors from output standard deviations.
#We must propagate the error because we fit the *inverse* parameters!
#The error = (error of the inverse)*(value**2)
for sd,value in zip(out.sd_beta, fit_out):
err_real=sd*(value**2)
fit_errors.append(err_real)
-
- def wlc_eval(x,params,pl_value,T):
+
+ def wlc_eval(x,params,pl_value,T):
'''
Evaluates the WLC function
'''
lambd, pii = params
else:
lambd = params
-
+
if pl_value:
pii=1/pl_value
-
+
Kb=(1.38065e-23) #boltzmann constant
therm=Kb*T #so we have thermal energy
-
+
return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
-
+
#STEP 3: plotting the fit
-
+
#obtain domain to plot the fit - from contact point to last_index plus 20 points
thule_index=last_index+10
if thule_index > len(xvector): #for rare cases in which we fit something at the END of whole curve.
xfit_chunk.reverse()
xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit_chunk]
xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
-
+
#the fitted curve: reflip, re-uncorrect
yfit=wlc_eval(xfit_chunk_corr_up, out.beta, pl_value,T)
yfit_down=[-y for y in yfit]
yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
-
+
if return_errors:
return fit_out, yfit_corr_down, xfit_chunk, fit_errors
else:
return fit_out, yfit_corr_down, xfit_chunk, None
-
+
def fjc_fit(self,clicked_points,xvector,yvector, pl_value, T=293, return_errors=False):
'''
Freely-jointed chain function
clicked_points[0] is the contact point (calculated or hand-clicked)
clicked_points[1] and [2] are edges of chunk
'''
-
+
#STEP 1: Prepare the vectors to apply the fit.
if pl_value is not None:
pl_value=pl_value/(10**9)
-
+
#indexes of the selected chunk
first_index=min(clicked_points[1].index, clicked_points[2].index)
last_index=max(clicked_points[1].index, clicked_points[2].index)
-
+
#getting the chunk and reverting it
xchunk,ychunk=xvector[first_index:last_index],yvector[first_index:last_index]
xchunk.reverse()
- ychunk.reverse()
+ ychunk.reverse()
#put contact point at zero and flip around the contact point (the fit wants a positive growth for extension and force)
xchunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xchunk]
ychunk_corr_up=[-(y-clicked_points[0].graph_coords[1]) for y in ychunk]
-
-
+
+
#make them arrays
xchunk_corr_up=scipy.array(xchunk_corr_up)
ychunk_corr_up=scipy.array(ychunk_corr_up)
-
-
+
+
#STEP 2: actually do the fit
-
+
#Find furthest point of chunk and add it a bit; the fit must converge
#from an excess!
xchunk_high=max(xchunk_corr_up)
xchunk_high+=(xchunk_high/10)
-
+
#Here are the linearized start parameters for the WLC.
#[lambd=1/Lo , pii=1/P]
-
+
p0=[(1/xchunk_high),(1/(3.5e-10))]
p0_plfix=[(1/xchunk_high)]
'''
hyperbolic cotangent
'''
return (np.exp(2*z)+1)/(np.exp(2*z)-1)
-
+
def x_fjc(params,f,T=T):
'''
fjc function for ODR fitting
lambd,pii=params
Kb=(1.38065e-23)
therm=Kb*T
-
+
#x=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
return x
-
+
def x_fjc_plfix(params,f,pl_value=pl_value,T=T):
'''
fjc function for ODR fitting
#y=(therm*pii/4.0) * (((1-(x*lambd))**-2) - 1 + (4*x*lambd))
x=(1/lambd)*(coth(f*(1/pii)/therm) - (therm*pii)/f)
return x
-
+
#make the ODR fit
realdata=scipy.odr.RealData(ychunk_corr_up,xchunk_corr_up)
if pl_value:
else:
model=scipy.odr.Model(x_fjc)
o = scipy.odr.ODR(realdata, model, p0)
-
+
o.set_job(fit_type=2)
out=o.run()
fit_out=[(1/i) for i in out.beta]
-
+
#Calculate fit errors from output standard deviations.
#We must propagate the error because we fit the *inverse* parameters!
#The error = (error of the inverse)*(value**2)
for sd,value in zip(out.sd_beta, fit_out):
err_real=sd*(value**2)
fit_errors.append(err_real)
-
- def fjc_eval(y,params,pl_value,T):
+
+ def fjc_eval(y,params,pl_value,T):
'''
Evaluates the WLC function
'''
lambd, pii = params
else:
lambd = params
-
+
if pl_value:
pii=1/pl_value
-
+
Kb=(1.38065e-23) #boltzmann constant
therm=Kb*T #so we have thermal energy
#return ( (therm*pii/4.0) * (((1-(x*lambd))**-2.0) - 1 + (4.0*x*lambd)) )
return (1/lambd)*(coth(y*(1/pii)/therm) - (therm*pii)/y)
-
-
+
+
#STEP 3: plotting the fit
#obtain domain to plot the fit - from contact point to last_index plus 20 points
thule_index=last_index+10
#or other buggy situations. Kludge to live with it now...
ychunk=yvector[:thule_index]
y_evalchunk=np.linspace(min(ychunk),max(ychunk),100)
-
+
yfit_down=[-y for y in y_evalchunk]
yfit_corr_down=[y+clicked_points[0].graph_coords[1] for y in yfit_down]
yfit_corr_down=scipy.array(yfit_corr_down)
-
+
#the fitted curve: reflip, re-uncorrect
xfit=fjc_eval(yfit_corr_down, out.beta, pl_value,T)
xfit=list(xfit)
xfit.reverse()
xfit_chunk_corr_up=[-(x-clicked_points[0].graph_coords[0]) for x in xfit]
-
+
#xfit_chunk_corr_up=scipy.array(xfit_chunk_corr_up)
#deltay=yfit_down[0]-yvector[clicked_points[0].index]
-
+
#This is a terrible, terrible kludge to find the point where it should normalize (and from where it should plot)
xxxdists=[]
for index in scipy.arange(1,len(xfit_chunk_corr_up),1):
- xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
+ xxxdists.append((clicked_points[0].graph_coords[0]-xfit_chunk_corr_up[index])**2)
normalize_index=xxxdists.index(min(xxxdists))
#End of kludge
-
+
deltay=yfit_down[normalize_index]-clicked_points[0].graph_coords[1]
yfit_corr_down=[y-deltay for y in yfit_down]
-
+
if return_errors:
return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], fit_errors
else:
return fit_out, yfit_corr_down[normalize_index+1:], xfit_chunk_corr_up[normalize_index+1:], None
-
-
+
+
def do_wlc(self,args):
'''
WLC
(fit.py plugin)
-
+
See the fit command
'''
self.do_fit(args)
-
+
def do_fjc(self,args):
'''
FJC
(fit.py plugin)
-
+
See the fit command
'''
self.do_fit(args)
-
+
def do_fit(self,args):
'''
FIT
First you have to click a contact point.
Then you have to click the two edges of the data you want to fit.
-
+
The fit function depends on the fit_function variable. You can set it with the command
"set fit_function wlc" or "set fit_function fjc" depending on the function you prefer.
-
- For WLC, the function is the simple polynomial worm-like chain as proposed by
- C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
+
+ For WLC, the function is the simple polynomial worm-like chain as proposed by
+ C.Bustamante, J.F.Marko, E.D.Siggia and S.Smith (Science. 1994
Sep 9;265(5178):1599-600.)
-
- For FJC, ref:
+
+ For FJC, ref:
C.Ray and B.B. Akhremitchev; http://www.chem.duke.edu/~boris/research/force_spectroscopy/fit_efjc.pdf
Arguments:
- pl=[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
+ pl=[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.
-
+ 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.
-
- noauto : allows for clicking the contact point by
+
+ noauto : allows for clicking the contact point by
hand (otherwise it is automatically estimated) the first time.
If subsequent measurements are made, the same contact point
clicked is used
-
+
reclick : redefines by hand the contact point, if noauto has been used before
but the user is unsatisfied of the previously choosen contact point.
---------
if ('t=' in arg[0:2]) or ('T=' in arg[0:2]):
t_expression=arg.split('=')
T=float(t_expression[1])
-
+
#use the currently displayed plot for the fit
displayed_plot=self._get_displayed_plot()
-
+
#handle contact point arguments correctly
if 'reclick' in args.split():
print 'Click contact point'
contact_point.absolute_coords=displayed_plot.vectors[1][0][cindex], displayed_plot.vectors[1][1][cindex]
contact_point.find_graph_coords(displayed_plot.vectors[1][0], displayed_plot.vectors[1][1])
contact_point.is_marker=True
-
+
print 'Click edges of chunk'
points=self._measure_N_points(N=2, whatset=1)
points=[contact_point]+points
print 'No recognized fit function defined!'
print 'Set your fit function to wlc or fjc.'
return
-
+
except:
print 'Fit not possible. Probably wrong interval -did you click two *different* points?'
return
-
+
#FIXME: print "Kuhn length" for FJC
print 'Fit function:',self.config['fit_function']
print 'Contour length: ',params[0]*(1.0e+9),' nm'
print name_of_charlength+': ',params[1]*(1.0e+9),' nm'
to_dump='persistent '+self.current.path+' '+str(params[1]*(1.0e+9))+' nm'
self.outlet.push(to_dump)
-
+
if fit_errors:
fit_nm=[i*(10**9) for i in fit_errors]
print 'Standard deviation (contour length)', fit_nm[0]
if len(fit_nm)>1:
print 'Standard deviation ('+name_of_charlength+')', fit_nm[1]
-
-
+
+
#add the clicked points in the final PlotObject
clickvector_x, clickvector_y=[], []
for item in points:
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
-
+
fitplot=copy.deepcopy(displayed_plot)
fitplot.add_set(xfit,yfit)
fitplot.add_set(clickvector_x,clickvector_y)
-
+
#FIXME: this colour/styles stuff must be solved at the root!
if fitplot.styles==[]:
fitplot.styles=[None,None,None,'scatter']
else:
fitplot.styles+=[None,'scatter']
-
+
if fitplot.colors==[]:
fitplot.colors=[None,None,None,None]
else:
fitplot.colors+=[None,None]
-
+
self._send_plot([fitplot])
-
-
+
+
#----------
-
-
+
+
def find_contact_point(self,plot=False):
'''
Finds the contact point on the curve.
-
+
The current algorithm (thanks to Francesco Musiani, francesco.musiani@unibo.it and Massimo Sandal) is:
- take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation
- fit the second half of the retraction curve to a line
- if the fit is not almost horizontal, take a smaller chunk and repeat
- otherwise, we have something horizontal
- so take the average of horizontal points and use it as a baseline
-
+
Then, start from the rise of the retraction curve and look at the first point below the
baseline.
-
+
FIXME: should be moved, probably to generalvclamp.py
'''
-
+
if not plot:
plot=self.plots[0]
-
+
outplot=self.subtract_curves(1)
xret=outplot.vectors[1][0]
ydiff=outplot.vectors[1][1]
yext=plot.vectors[0][1]
xret2=plot.vectors[1][0]
yret=plot.vectors[1][1]
-
+
#taking care of the picoforce trigger bug: we exclude portions of the curve that have too much
#standard deviation. yes, a lot of magic is here.
monster=True
else: #move away from the monster
monlength-=int(len(xret)/50)
finalength-=int(len(xret)/50)
-
-
+
+
#take half of the thing
endlength=int(len(xret)/2)
-
+
ok=False
-
+
while not ok:
xchunk=yext[endlength:monlength]
ychunk=yext[endlength:monlength]
if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) :
endlength+=10
else:
- ok=True
-
-
+ ok=True
+
+
ymean=np.mean(ychunk) #baseline
-
+
index=0
point = ymean+1
-
+
#find the first point below the calculated baseline
while point > ymean:
try:
point=yret[index]
- index+=1
+ index+=1
except IndexError:
#The algorithm didn't find anything below the baseline! It should NEVER happen
- index=0
+ index=0
return index
-
+
return index
-
-
-
+
+
+
def find_contact_point2(self, debug=False):
'''
TO BE DEVELOPED IN THE FUTURE
Finds the contact point on the curve.
-
+
FIXME: should be moved, probably to generalvclamp.py
'''
-
+
#raw_plot=self.current.curve.default_plots()[0]
raw_plot=self.plots[0]
'''xext=self.plots[0].vectors[0][0]
yext=raw_plot.vectors[0][1]
xret2=raw_plot.vectors[1][0]
yret=raw_plot.vectors[1][1]
-
+
first_point=[xext[0], yext[0]]
last_point=[xext[-1], yext[-1]]
-
+
#regr=scipy.polyfit(first_point, last_point,1)[0:2]
diffx=abs(first_point[0]-last_point[0])
diffy=abs(first_point[1]-last_point[1])
-
+
#using polyfit results in numerical errors. good old algebra.
a=diffy/diffx
b=first_point[1]-(a*first_point[0])
baseline=scipy.polyval((a,b), xext)
-
+
ysub=[item-basitem for item,basitem in zip(yext,baseline)]
-
+
contact=ysub.index(min(ysub))
-
+
return xext,ysub,contact
-
+
#now, exploit a ClickedPoint instance to calculate index...
dummy=ClickedPoint()
dummy.absolute_coords=(x_intercept,y_intercept)
dummy.find_graph_coords(xret2,yret)
-
+
if debug:
return dummy.index, regr, regr_contact
else:
return dummy.index
-
-
+
+
def x_do_contact(self,args):
'''
DEBUG COMMAND to be activated in the future
'''
xext,ysub,contact=self.find_contact_point2(debug=True)
-
+
contact_plot=self.plots[0]
contact_plot.add_set(xext,ysub)
contact_plot.add_set([xext[contact]],[self.plots[0].vectors[0][1][contact]])
contact_plot.styles=[None,None,None,'scatter']
self._send_plot([contact_plot])
return
-
-
+
+
index,regr,regr_contact=self.find_contact_point2(debug=True)
print regr
print regr_contact
#nc_line=[(item*regr[0])+regr[1] for item in x_nc]
nc_line=scipy.polyval(regr,xret)
c_line=scipy.polyval(regr_contact,xret)
-
-
+
+
contact_plot=self.current.curve.default_plots()[0]
contact_plot.add_set(xret, nc_line)
contact_plot.add_set(xret, c_line)
#contact_plot.styles.append(None)
contact_plot.destination=1
self._send_plot([contact_plot])
-
\ No newline at end of file
-#!/usr/bin/env python
-
'''
FLATFILTS
class flatfiltsCommands:
-
+
def _plug_init(self):
#configurate convfilt variables
convfilt_configurator=ConvfiltConfig()
-
+
#different OSes have different path conventions
if self.config['hookedir'][0]=='/':
slash='/' #a Unix or Unix-like system
else:
slash='\\' #it's a drive letter, we assume it's Windows
-
+
self.convfilt_config=convfilt_configurator.load_config(self.config['hookedir']+slash+'convfilt.conf')
-
+
def do_flatfilt(self,args):
'''
FLATFILT
median_filter=7
min_npks=4
min_deviation=9
-
+
args=args.split(' ')
if len(args) == 2:
min_npks=int(args[0])
min_deviation=int(args[1])
else:
pass
-
+
print 'Processing playlist...'
notflat_list=[]
-
+
c=0
for item in self.current_list:
c+=1
-
+
try:
notflat=self.has_features(item, median_filter, min_npks, min_deviation)
print 'Curve',item.path, 'is',c,'of',len(self.current_list),': features are ',notflat
except:
notflat=False
print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
-
+
if notflat:
item.features=notflat
item.curve=None #empty the item object, to further avoid memory leak
notflat_list.append(item)
-
+
if len(notflat_list)==0:
print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
return
self.current_list=notflat_list
self.current=self.current_list[self.pointer]
self.do_plot(0)
-
+
def has_features(self,item,median_filter,min_npks,min_deviation):
'''
decides if a curve is flat enough to be rejected from analysis: it sees if there
are at least min_npks points that are higher than min_deviation times the absolute value
of noise.
-
+
Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
'''
retvalue=False
-
- item.identify(self.drivers)
+
+ item.identify(self.drivers)
#we assume the first is the plot with the force curve
#do the median to better resolve features from noise
flat_plot=self.plotmanip_median(item.curve.default_plots()[0], item, customvalue=median_filter)
- flat_vects=flat_plot.vectors
+ flat_vects=flat_plot.vectors
item.curve.close_all()
#needed to avoid *big* memory leaks!
del item.curve
del item
-
- #absolute value of derivate
+
+ #absolute value of derivate
yretdiff=diff(flat_vects[1][1])
yretdiff=[abs(value) for value in yretdiff]
#average of derivate values
c_pks+=1
else:
break
-
+
if c_pks>=min_npks:
retvalue = c_pks
-
+
del flat_plot, flat_vects, yretdiff
-
+
return retvalue
################################################################
#-----CONVFILT-------------------------------------------------
#-----Convolution-based peak recognition and filtering.
#Requires the libpeakspot.py library
-
+
def has_peaks(self, plot, abs_devs=None):
'''
Finds peak position in a force curve.
'''
if abs_devs==None:
abs_devs=self.convfilt_config['mindeviation']
-
-
+
+
xret=plot.vectors[1][0]
yret=plot.vectors[1][1]
#Calculate convolution.
convoluted=lps.conv_dx(yret, self.convfilt_config['convolution'])
-
+
#surely cut everything before the contact point
cut_index=self.find_contact_point(plot)
#cut even more, before the blind window
blind_index+=1
cut_index+=blind_index
#do the dirty convolution-peak finding stuff
- noise_level=lps.noise_absdev(convoluted[cut_index:], self.convfilt_config['positive'], self.convfilt_config['maxcut'], self.convfilt_config['stable'])
- above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs)
+ noise_level=lps.noise_absdev(convoluted[cut_index:], self.convfilt_config['positive'], self.convfilt_config['maxcut'], self.convfilt_config['stable'])
+ above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs)
peak_location,peak_size=lps.find_peaks(above,seedouble=self.convfilt_config['seedouble'])
-
+
#take the maximum
for i in range(len(peak_location)):
peak=peak_location[i]
maxpk=min(yret[peak-10:peak+10])
index_maxpk=yret[peak-10:peak+10].index(maxpk)+(peak-10)
peak_location[i]=index_maxpk
-
+
return peak_location,peak_size
-
-
+
+
def exec_has_peaks(self,item,abs_devs):
'''
encapsulates has_peaks for the purpose of correctly treating the curve objects in the convfilt loop,
to avoid memory leaks
'''
- item.identify(self.drivers)
+ item.identify(self.drivers)
#we assume the first is the plot with the force curve
plot=item.curve.default_plots()[0]
-
+
if 'flatten' in self.config['plotmanips']:
#If flatten is present, use it for better recognition of peaks...
flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
plot=flatten(plot, item, customvalue=1)
-
+
peak_location,peak_size=self.has_peaks(plot,abs_devs)
#close all open files
item.curve.close_all()
del item.curve
del item
return peak_location, peak_size
-
+
#------------------------
#------commands----------
- #------------------------
+ #------------------------
def do_peaks(self,args):
'''
PEAKS
'''
if len(args)==0:
args=self.convfilt_config['mindeviation']
-
+
try:
abs_devs=float(args)
except:
print 'Wrong argument, using config value'
abs_devs=float(self.convfilt_config['mindeviation'])
-
+
defplots=self.current.curve.default_plots()[0] #we need the raw, uncorrected plots
-
+
if 'flatten' in self.config['plotmanips']:
flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
defplots=flatten(defplots, self.current)
else:
print 'You have the flatten plot manipulator not loaded. Enabling it could give you better results.'
-
+
peak_location,peak_size=self.has_peaks(defplots,abs_devs)
print 'Found '+str(len(peak_location))+' peaks.'
to_dump='peaks '+self.current.path+' '+str(len(peak_location))
self.outlet.push(to_dump)
#print peak_location
-
+
#if no peaks, we have nothing to plot. exit.
if len(peak_location)==0:
return
-
+
#otherwise, we plot the peak locations.
xplotted_ret=self.plots[0].vectors[1][0]
yplotted_ret=self.plots[0].vectors[1][1]
xgood=[xplotted_ret[index] for index in peak_location]
ygood=[yplotted_ret[index] for index in peak_location]
-
+
recplot=self._get_displayed_plot()
recplot.vectors.append([xgood,ygood])
if recplot.styles==[]:
else:
recplot.styles+=['scatter']
recplot.colors+=[None]
-
+
self._send_plot([recplot])
-
+
def do_convfilt(self,args):
'''
CONVFILT
If called without arguments, it uses default values.
'''
-
+
min_npks=self.convfilt_config['minpeaks']
min_deviation=self.convfilt_config['mindeviation']
-
+
args=args.split(' ')
if len(args) == 2:
min_npks=int(args[0])
min_deviation=int(args[1])
else:
pass
-
+
print 'Processing playlist...'
print '(Please wait)'
notflat_list=[]
-
+
c=0
for item in self.current_list:
c+=1
-
- try:
+
+ try:
peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
if len(peak_location)>=min_npks:
isok='+'
except:
peak_location,peak_size=[],[]
print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
-
+
if len(peak_location)>=min_npks:
item.peak_location=peak_location
item.peak_size=peak_size
item.curve=None #empty the item object, to further avoid memory leak
notflat_list.append(item)
-
+
#Warn that no flattening had been done.
if not ('flatten' in self.config['plotmanips']):
print 'Flatten manipulator was not found. Processing was done without flattening.'
print 'Try to enable it in your configuration file for better results.'
-
+
if len(notflat_list)==0:
print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
return
self.current_list=notflat_list
self.current=self.current_list[self.pointer]
self.do_plot(0)
-
-
+
+
def do_setconv(self,args):
'''
SETCONV
print 'This is not an internal convfilt variable!'
print 'Run "setconv" without arguments to see a list of defined variables.'
return
-
+
if len(args)==1:
print self.convfilt_config[args[0]]
elif len(args)>1:
class ConvfiltConfig:
'''
Handling of convfilt configuration file
-
+
Mostly based on the simple-yet-useful examples of the Python Library Reference
about xml.dom.minidom
-
+
FIXME: starting to look a mess, should require refactoring
'''
-
+
def __init__(self):
self.config={}
-
-
+
+
def load_config(self, filename):
- myconfig=file(filename)
+ myconfig=file(filename)
#the following 3 lines are needed to strip newlines. otherwise, since newlines
#are XML elements too, the parser would read them (and re-save them, multiplying
#newlines...)
the_file=myconfig.read()
the_file_lines=the_file.split('\n')
the_file=''.join(the_file_lines)
-
- self.config_tree=xml.dom.minidom.parseString(the_file)
-
+
+ self.config_tree=xml.dom.minidom.parseString(the_file)
+
def getText(nodelist):
#take the text from a nodelist
#from Python Library Reference 13.7.2
if node.nodeType == node.TEXT_NODE:
rc += node.data
return rc
-
+
def handleConfig(config):
noiseabsdev_elements=config.getElementsByTagName("noise_absdev")
convfilt_elements=config.getElementsByTagName("convfilt")
handleAbsdev(noiseabsdev_elements)
handleConvfilt(convfilt_elements)
-
+
def handleAbsdev(noiseabsdev_elements):
for element in noiseabsdev_elements:
for attribute in element.attributes.keys():
self.config[attribute]=element.getAttribute(attribute)
-
+
def handleConvfilt(convfilt_elements):
for element in convfilt_elements:
for attribute in element.attributes.keys():
self.config[attribute]=element.getAttribute(attribute)
-
+
handleConfig(self.config_tree)
#making items in the dictionary machine-readable
for item in self.config.keys():
self.config[item]=eval(self.config[item])
except NameError: #if it's an unreadable string, keep it as a string
pass
-
- return self.config
\ No newline at end of file
+
+ return self.config
-#!/usr/bin/env python
-
'''
GENERALCLAMP.py
Plugin regarding general force clamp measurements
'''
from libhooke import WX_GOOD, ClickedPoint
-import wxversion\r
+import wxversion
import libhookecurve as lhc
wxversion.select(WX_GOOD)
-from wx import PostEvent\r
+from wx import PostEvent
+
+class generalclampCommands:
-class generalclampCommands:\r
-\r
def plotmanip_clamp(self, plot, current, customvalue=False):
'''
- Handles some viewing options for the "force clamp" data format, depending on the state of these configuration variables:\r
- (1) If self.config['fc_showphase'] != 0, the 'phase' data column (i.e. the 2nd) is shown in the 0th graph (else it isn't)\r
- (2) If self.config['fc_showimposed'] != 0, the 'imposed deflection' data column (i.e. the 5th) is shown in the 1st graph (else it isn't)\r
+ Handles some viewing options for the "force clamp" data format, depending on the state of these configuration variables:
+ (1) If self.config['fc_showphase'] != 0, the 'phase' data column (i.e. the 2nd) is shown in the 0th graph (else it isn't)
+ (2) If self.config['fc_showimposed'] != 0, the 'imposed deflection' data column (i.e. the 5th) is shown in the 1st graph (else it isn't)
(3) If self.config['fc_interesting'] == 0, the entire curve is shown in the graphs; if it has a non-zero value N, only phase N is shown.
-\r
- NOTE - my implementation of point(3) feels quite awkward - someone smarter than me plz polish that!\r
-\r
+
+ NOTE - my implementation of point(3) feels quite awkward - someone smarter than me plz polish that!
+
'''
-
+
#not a fclamp curve...
if current.curve.experiment != 'clamp':
- return plot\r
-\r
- if self.config['fc_interesting'] != 0 and plot.destination==0:\r
- lower = int((self.config['fc_interesting'])-1)\r
- upper = int((self.config['fc_interesting'])+1)\r
- trim = current.curve.trimindexes()[lower:upper]\r
- newtime = []\r
- newzpiezo = []\r
- newphase = []\r
- for x in range(trim[0],trim[1]):\r
- newtime.append(self.plots[0].vectors[0][0][x])\r
- newzpiezo.append(self.plots[0].vectors[0][1][x])\r
- newphase.append(self.plots[0].vectors[1][1][x])\r
- self.plots[0].vectors[0][0] = newtime\r
- self.plots[0].vectors[0][1] = newzpiezo\r
- self.plots[0].vectors[1][0] = newtime\r
- self.plots[0].vectors[1][1] = newphase\r
-\r
- if self.config['fc_interesting'] != 0 and plot.destination==1:\r
- lower = int((self.config['fc_interesting'])-1)\r
- upper = int((self.config['fc_interesting'])+1)\r
- trim = current.curve.trimindexes()[lower:upper]\r
- newtime = []\r
- newdefl = []\r
- newimposed = []\r
- for x in range(trim[0],trim[1]):\r
- newtime.append(self.plots[1].vectors[0][0][x])\r
- newdefl.append(self.plots[1].vectors[0][1][x])\r
- newimposed.append(self.plots[1].vectors[1][1][x])\r
- self.plots[1].vectors[0][0] = newtime\r
- self.plots[1].vectors[0][1] = newdefl\r
- self.plots[1].vectors[1][0] = newtime\r
- self.plots[1].vectors[1][1] = newimposed \r
- \r
- if self.config['fc_showphase'] == 0 and plot.destination==0:\r
- self.plots[0].remove_set(1)\r
- \r
- if self.config['fc_showimposed'] == 0 and plot.destination==1:\r
- self.plots[1].remove_set(1)\r
- \r
+ return plot
+
+ if self.config['fc_interesting'] != 0 and plot.destination==0:
+ lower = int((self.config['fc_interesting'])-1)
+ upper = int((self.config['fc_interesting'])+1)
+ trim = current.curve.trimindexes()[lower:upper]
+ newtime = []
+ newzpiezo = []
+ newphase = []
+ for x in range(trim[0],trim[1]):
+ newtime.append(self.plots[0].vectors[0][0][x])
+ newzpiezo.append(self.plots[0].vectors[0][1][x])
+ newphase.append(self.plots[0].vectors[1][1][x])
+ self.plots[0].vectors[0][0] = newtime
+ self.plots[0].vectors[0][1] = newzpiezo
+ self.plots[0].vectors[1][0] = newtime
+ self.plots[0].vectors[1][1] = newphase
+
+ if self.config['fc_interesting'] != 0 and plot.destination==1:
+ lower = int((self.config['fc_interesting'])-1)
+ upper = int((self.config['fc_interesting'])+1)
+ trim = current.curve.trimindexes()[lower:upper]
+ newtime = []
+ newdefl = []
+ newimposed = []
+ for x in range(trim[0],trim[1]):
+ newtime.append(self.plots[1].vectors[0][0][x])
+ newdefl.append(self.plots[1].vectors[0][1][x])
+ newimposed.append(self.plots[1].vectors[1][1][x])
+ self.plots[1].vectors[0][0] = newtime
+ self.plots[1].vectors[0][1] = newdefl
+ self.plots[1].vectors[1][0] = newtime
+ self.plots[1].vectors[1][1] = newimposed
+
+ if self.config['fc_showphase'] == 0 and plot.destination==0:
+ self.plots[0].remove_set(1)
+
+ if self.config['fc_showimposed'] == 0 and plot.destination==1:
+ self.plots[1].remove_set(1)
+
return plot
-
+
def do_time(self,args):
'''
Measures the time difference (in seconds) between two points
----
Syntax: time
'''
- if self.current.curve.experiment == 'clamp':\r
- time=self._delta(set=0)[0]\r
+ if self.current.curve.experiment == 'clamp':
+ time=self._delta(set=0)[0]
print str(time*1000)+' ms'
else:
print 'This command makes no sense for a non-force clamp experiment.'
-
+
def do_zpiezo(self,args):
'''
Measures the zpiezo difference (in nm) between two points
Syntax: zpiezo
'''
if self.current.curve.experiment == 'clamp':
- zpiezo=self._delta(set=0)[2]\r
+ zpiezo=self._delta(set=0)[2]
print str(zpiezo*(10**9))+' nm'
else:
print 'This command makes no sense for a non-force clamp experiment.'
-
+
def do_defl(self,args):
'''
Measures the deflection difference (in nm) between two points
print str(defl*(10**12))+' pN'
else:
print 'This command makes no sense for a non-force clamp experiment.'
-
+
def do_step(self,args):
'''
Measures the length and time duration of a time-Z step
dt=abs(points[1].graph_coords[0]-points[0].graph_coords[0])
print 'dZ: ',dz,' nm'
print 'dT: ',dt,' s'
-
+
+ else:
+ print 'This command makes no sense for a non-force clamp experiment.'
+
+ def do_fcfilt(self,args):
+ '''
+ Filters out featureless force clamp curves of the current playlist.
+ It's very similar to 'flatfilt' for velocity clamp curves.
+ Creates a new playlist only containing non-empty curves.
+
+ WARNING - Only works if you set an appropriate fc_interesting config variable!
+ WARNING - arguments are NOT optional at the moment!
+
+ Syntax: fcfilt maxretraction(nm) mindeviation (pN)
+
+ Suggested values for an (i27)8 experiment with our setup are 200nm and 10-15 pN
+ '''
+
+ if self.config['fc_interesting'] == 0:
+ print 'You must specify the phase of interest (using set fc_interesing X) prior to running fcfilt!'
+ return
+
+ maxretraction=0
+ threshold=0
+ args=args.split(' ')
+ if len(args)==2:
+ maxretraction=int(args[0])
+ threshold=int(args[1])
+ else:
+ print 'Arguments are not optional for fcfilt. You should pass two numbers:'
+ print '(1) the maximum plausible piezo retraction in NANOMETERS (e.g. the length of the protein)'
+ print "(2) the threshold, in PICONEWTONS. If signal deviates from imposed more than this, it's an event"
+ return
+
+
+ print 'Processing playlist... go get yourself a cup of coffee.'
+ notflat_list=[]
+
+ c=0
+
+ for item in self.current_list:
+ c+=1
+ try:
+ notflat=self.has_stuff(item,maxretraction,threshold)
+ print 'Curve',item.path,'is',c,'of',len(self.current_list),'--->Has Stuff =',notflat
+ except:
+ notflat=False
+ print 'Curve',item.path,'is',c,'of',len(self.current_list),'--->could not be processed'
+ if notflat:
+ item.features=notflat
+ item.curve=None
+ notflat_list.append(item)
+
+ if len(notflat_list)==0:
+ print 'Nothing interesting here. Reconsider either your filtering criteria or your experimental data'
+ return
+ else:
+ print 'Found',len(notflat_list),'potentially interesting curves.'
+ print 'Regenerating Playlist...'
+ self.pointer=0
+ self.current_list=notflat_list
+ self.current=self.current_list[self.pointer]
+ self.do_plot(0)
+
+ def has_stuff(self,item,maxretraction,threshold):
+ '''
+ Decides whether a curve has some features in the interesting phase.
+ Algorithm:
+ - clip the interesting phase portion of the curve.
+ - discard the first 20 milliseconds (this is due to a quirk of our hardware).
+ - look at the zpiezo plot and note down when (if) retratcs more than [maxretraction] nm away from the first point.
+ - clip off any data after this point, with an excess of 100 points (again, an hardware quirk)
+ - if the remainder is less than 100 points, ditch the curve.
+ - now look at the deflection plot and check if there are points more than [threshold] pN over the 'flat zone'.
+ - if you find such points, bingo!
+ '''
+
+ item.identify(self.drivers)
+
+ lower = int((self.config['fc_interesting'])-1)
+ upper = int((self.config['fc_interesting'])+1)
+ trim_idxs = item.curve.trimindexes()[lower:upper]
+ lo=trim_idxs[0]+20 #clipping the first 20 points off...
+ hi=trim_idxs[1]
+ trimmed_zpiezo=item.curve.default_plots()[0].vectors[0][1][lo:hi]
+ trimmed_defl=item.curve.default_plots()[1].vectors[0][1][lo:hi]
+ trimmed_imposed=item.curve.default_plots()[1].vectors[1][1][lo:hi]
+ imposed=trimmed_imposed[21] #just to match the 20-pts clipping...
+
+ item.curve.close_all()
+ del item.curve
+ del item
+
+ starting_z=trimmed_zpiezo[0]
+ plausible=starting_z-(maxretraction*1e-9)
+ det_trim=0
+ while trimmed_zpiezo[det_trim]>plausible:
+ det_trim+=1
+ if det_trim >= len(trimmed_zpiezo): #breaking cycles makes me shiver...
+ det_trim=len(trimmed_zpiezo) #but I cannot think of anything better now.
+ break
+ further_trim=det_trim-100
+ if further_trim<100:
+ return False
+ trimmed_defl=trimmed_defl[:further_trim]
+
+ trimmed_defl.sort()
+ ninetypercent=int(0.9*len(trimmed_defl))
+ j=0
+ sum=0
+ for j in trimmed_defl[:ninetypercent]:
+ sum+=j
+ avg=float(sum/ninetypercent)
+ sweetspot=float(avg+(threshold*1e-12))
+ if trimmed_defl[-1]>sweetspot:
+ flag=True
else:
- print 'This command makes no sense for a non-force clamp experiment.'\r
-\r
- def do_fcfilt(self,args):\r
- '''\r
- Filters out featureless force clamp curves of the current playlist.\r
- It's very similar to 'flatfilt' for velocity clamp curves.\r
- Creates a new playlist only containing non-empty curves.\r
-\r
- WARNING - Only works if you set an appropriate fc_interesting config variable!\r
- WARNING - arguments are NOT optional at the moment!\r
-\r
- Syntax: fcfilt maxretraction(nm) mindeviation (pN)\r
-\r
- Suggested values for an (i27)8 experiment with our setup are 200nm and 10-15 pN\r
- '''\r
-\r
- if self.config['fc_interesting'] == 0:\r
- print 'You must specify the phase of interest (using set fc_interesing X) prior to running fcfilt!'\r
- return\r
- \r
- maxretraction=0\r
- threshold=0\r
- args=args.split(' ')\r
- if len(args)==2:\r
- maxretraction=int(args[0])\r
- threshold=int(args[1])\r
- else:\r
- print 'Arguments are not optional for fcfilt. You should pass two numbers:'\r
- print '(1) the maximum plausible piezo retraction in NANOMETERS (e.g. the length of the protein)'\r
- print "(2) the threshold, in PICONEWTONS. If signal deviates from imposed more than this, it's an event"\r
- return\r
- \r
-\r
- print 'Processing playlist... go get yourself a cup of coffee.'\r
- notflat_list=[]\r
-\r
- c=0\r
-\r
- for item in self.current_list:\r
- c+=1\r
- try:\r
- notflat=self.has_stuff(item,maxretraction,threshold)\r
- print 'Curve',item.path,'is',c,'of',len(self.current_list),'--->Has Stuff =',notflat\r
- except:\r
- notflat=False\r
- print 'Curve',item.path,'is',c,'of',len(self.current_list),'--->could not be processed'\r
- if notflat:\r
- item.features=notflat\r
- item.curve=None\r
- notflat_list.append(item)\r
-\r
- if len(notflat_list)==0:\r
- print 'Nothing interesting here. Reconsider either your filtering criteria or your experimental data'\r
- return\r
- else:\r
- print 'Found',len(notflat_list),'potentially interesting curves.'\r
- print 'Regenerating Playlist...'\r
- self.pointer=0\r
- self.current_list=notflat_list\r
- self.current=self.current_list[self.pointer]\r
- self.do_plot(0)\r
-\r
- def has_stuff(self,item,maxretraction,threshold):\r
- '''\r
- Decides whether a curve has some features in the interesting phase.\r
- Algorithm:\r
- - clip the interesting phase portion of the curve.\r
- - discard the first 20 milliseconds (this is due to a quirk of our hardware).\r
- - look at the zpiezo plot and note down when (if) retratcs more than [maxretraction] nm away from the first point.\r
- - clip off any data after this point, with an excess of 100 points (again, an hardware quirk)\r
- - if the remainder is less than 100 points, ditch the curve.\r
- - now look at the deflection plot and check if there are points more than [threshold] pN over the 'flat zone'.\r
- - if you find such points, bingo! \r
- '''\r
-\r
- item.identify(self.drivers)\r
- \r
- lower = int((self.config['fc_interesting'])-1)\r
- upper = int((self.config['fc_interesting'])+1)\r
- trim_idxs = item.curve.trimindexes()[lower:upper]\r
- lo=trim_idxs[0]+20 #clipping the first 20 points off...\r
- hi=trim_idxs[1]\r
- trimmed_zpiezo=item.curve.default_plots()[0].vectors[0][1][lo:hi]\r
- trimmed_defl=item.curve.default_plots()[1].vectors[0][1][lo:hi]\r
- trimmed_imposed=item.curve.default_plots()[1].vectors[1][1][lo:hi]\r
- imposed=trimmed_imposed[21] #just to match the 20-pts clipping...\r
- \r
- item.curve.close_all()\r
- del item.curve\r
- del item\r
-\r
- starting_z=trimmed_zpiezo[0]\r
- plausible=starting_z-(maxretraction*1e-9)\r
- det_trim=0\r
- while trimmed_zpiezo[det_trim]>plausible:\r
- det_trim+=1\r
- if det_trim >= len(trimmed_zpiezo): #breaking cycles makes me shiver...\r
- det_trim=len(trimmed_zpiezo) #but I cannot think of anything better now.\r
- break\r
- further_trim=det_trim-100\r
- if further_trim<100:\r
- return False\r
- trimmed_defl=trimmed_defl[:further_trim]\r
-\r
- trimmed_defl.sort()\r
- ninetypercent=int(0.9*len(trimmed_defl))\r
- j=0\r
- sum=0\r
- for j in trimmed_defl[:ninetypercent]:\r
- sum+=j\r
- avg=float(sum/ninetypercent)\r
- sweetspot=float(avg+(threshold*1e-12))\r
- if trimmed_defl[-1]>sweetspot:\r
- flag=True\r
- else:\r
- flag=False\r
-\r
- del trimmed_defl,trimmed_zpiezo,trimmed_imposed \r
-\r
- return flag \r
-
\ No newline at end of file
+ flag=False
+
+ del trimmed_defl,trimmed_zpiezo,trimmed_imposed
+
+ return flag
+
-#!/usr/bin/env python
-
'''
generaltccd.py
'''
class generaltccdCommands:
-
+
def plotmanip_threshold(self, plot, current, customvalue=False):
'''
Cuts from the plot everything below the threshold.
Set the threshold with "set tccd_threshold"
'''
-
+
if current.curve.experiment != 'smfluo':
return plot
-
+
if not self.config['tccd_threshold'] and (not customvalue):
return plot
-
+
if customvalue:
thresh=customvalue
else:
thresh=self.config['tccd_threshold']
-
+
for set in plot.vectors:
newy=[]
for value in set[1]:
newy.append(0)
else:
newy.append(value)
-
+
set[1]=newy
-
+
return plot
-
+
def plotmanip_coincident(self,plot,current, customvalue=False):
'''
'''
if current.curve.experiment != 'smfluo':
return plot
-
+
if not self.config['tccd_coincident'] and (not customvalue):
return plot
-
+
newred=[]
newblue=[]
for index in range(len(plot.vectors[0][1])):
else:
newred.append(0)
newblue.append(0)
-
+
plot.vectors[0][1]=newred
plot.vectors[1][1]=newblue
-
- return plot
\ No newline at end of file
+
+ return plot
-#!/usr/bin/env python
-
'''
generalvclamp.py
class generalvclampCommands:
-
+
def _plug_init(self):
self.basecurrent=None
self.basepoints=None
self.autofile=''
-
+
def do_distance(self,args):
'''
DISTANCE
Measure the force difference (in pN) between two points
---------------
Syntax: force
- '''
+ '''
if self.current.curve.experiment == 'clamp':
print 'This command makes no sense for a force clamp experiment.'
return
print str(dy*(10**12))+' pN'
to_dump='force '+self.current.path+' '+str(dy*(10**12))+' pN'
self.outlet.push(to_dump)
-
-
+
+
def do_forcebase(self,args):
'''
FORCEBASE
(generalvclamp.py)
Measures the difference in force (in pN) between a point and a baseline
took as the average between two points.
-
+
The baseline is fixed once for a given curve and different force measurements,
unless the user wants it to be recalculated
------------
'''
rebase=False #if true=we select rebase
maxpoint=False #if true=we measure the maximum peak
-
+
plot=self._get_displayed_plot()
whatset=1 #fixme: for all sets
if 'rebase' in args or (self.basecurrent != self.current.path):
rebase=True
if 'max' in args:
maxpoint=True
-
+
if rebase:
print 'Select baseline'
self.basepoints=self._measure_N_points(N=2, whatset=whatset)
self.basecurrent=self.current.path
-
+
if maxpoint:
print 'Select two points'
points=self._measure_N_points(N=2, whatset=whatset)
points=self._measure_N_points(N=1, whatset=whatset)
#whatplot=points[0].dest
y=points[0].graph_coords[1]
-
+
#fixme: code duplication
boundaries=[self.basepoints[0].index, self.basepoints[1].index]
boundaries.sort()
to_average=plot.vectors[whatset][1][boundaries[0]:boundaries[1]] #y points to average
-
+
avg=np.mean(to_average)
forcebase=abs(y-avg)
print str(forcebase*(10**12))+' pN'
Multiplies all the Y values of an SMFS curve by a value stored in the 'force_multiplier'\r
configuration variable. Useful for calibrations and other stuff.
'''
-
+
#not a smfs curve...
if current.curve.experiment != 'smfs':
return plot
-
+
#only one set is present...
if len(self.plots[0].vectors) != 2:
return plot
-
+
#multiplier is 1...
if (self.config['force_multiplier']==1):
return plot\r
plot.vectors[1][1][i]=plot.vectors[1][1][i]*self.config['force_multiplier']\r
\r
return plot \r
-
-
+
+
def plotmanip_flatten(self, plot, current, customvalue=False):
'''
Subtracts a polynomial fit to the non-contact part of the curve, as to flatten it.
- the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
+ the best polynomial fit is chosen among polynomials of degree 1 to n, where n is
given by the configuration file or by the customvalue.
-
+
customvalue= int (>0) --> starts the function even if config says no (default=False)
'''
-
+
#not a smfs curve...
if current.curve.experiment != 'smfs':
return plot
-
+
#only one set is present...
if len(self.plots[0].vectors) != 2:
return plot
-
+
#config is not flatten, and customvalue flag is false too
if (not self.config['flatten']) and (not customvalue):
return plot
-
+
max_exponent=12
delta_contact=0
-
+
if customvalue:
max_cycles=customvalue
else:
max_cycles=self.config['flatten'] #Using > 1 usually doesn't help and can give artefacts. However, it could be useful too.
-
+
contact_index=self.find_contact_point()
-
+
valn=[[] for item in range(max_exponent)]
yrn=[0.0 for item in range(max_exponent)]
errn=[0.0 for item in range(max_exponent)]
-
+
for i in range(int(max_cycles)):
-
+
x_ext=plot.vectors[0][0][contact_index+delta_contact:]
y_ext=plot.vectors[0][1][contact_index+delta_contact:]
x_ret=plot.vectors[1][0][contact_index+delta_contact:]
return plot
best_exponent=errn.index(min(errn))
-
+
#extension
ycorr_ext=y_ext-yrn[best_exponent]+y_ext[0] #noncontact part
- yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
+ yjoin_ext=np.array(plot.vectors[0][1][0:contact_index+delta_contact]) #contact part
#retraction
ycorr_ret=y_ret-yrn[best_exponent]+y_ext[0] #noncontact part
yjoin_ret=np.array(plot.vectors[1][1][0:contact_index+delta_contact]) #contact part
-
+
ycorr_ext=np.concatenate((yjoin_ext, ycorr_ext))
ycorr_ret=np.concatenate((yjoin_ret, ycorr_ret))
-
+
plot.vectors[0][1]=list(ycorr_ext)
plot.vectors[1][1]=list(ycorr_ret)
-
+
return plot
-
+
#---SLOPE---
def do_slope(self,args):
'''
clickedpoints=[]
points=self._measure_N_points(N=1,whatset=1)
clickedpoints=[points[0].index-fitspan,points[0].index]
-
+
# Calls the function linefit_between
parameters=[0,0,[],[]]
try:
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=[]
x=0
for x in xtoplot:
ytoplot.append((x*parameters[0])+parameters[1])
-
+
clickvector_x, clickvector_y=[], []
for item in points:
clickvector_x.append(item.graph_coords[0])
clickvector_y.append(item.graph_coords[1])
lineplot=self._get_displayed_plot(0) #get topmost displayed plot
-
+
lineplot.add_set(xtoplot,ytoplot)
lineplot.add_set(clickvector_x, clickvector_y)
-
-
+
+
if lineplot.styles==[]:
lineplot.styles=[None,None,None,'scatter']
else:
lineplot.styles=[None,None,None,None]
else:
lineplot.colors+=[None,None]
-
-
+
+
self._send_plot([lineplot])
def linefit_between(self,index1,index2,whatset=1):
# Translates the indexes into two vectors containing the x,y data to fit
xtofit=self.plots[0].vectors[whatset][0][index1:index2]
ytofit=self.plots[0].vectors[whatset][1][index1:index2]
-
+
# Does the actual linear fitting (simple least squares with numpy.polyfit)
linefit=[]
linefit=np.polyfit(xtofit,ytofit,1)
return (linefit[0],linefit[1],xtofit,ytofit)
-
-
-
-#!/usr/bin/env python
-
'''
COMMAND MACRO PLUGIN FOR HOOKE
if args=='stop':
self.pause=0
self.prompt=self.auxprompt
- if len(self.currentmacro) != 0:
+ if len(self.currentmacro) != 0:
answer=linput.safeinput('Do you want to save this macro? ',['y'])
if answer[0].lower() == 'y':
self.do_savemacro('')
self.pause=0
self.collect()
else:
- if len(self.currentmacro) != 0:
+ if len(self.currentmacro) != 0:
answer=linput.safeinput('Another macro is already beign recorded\nDo you want to save it?',['y'])
if answer[0].lower() == 'y':
self.do_savemacro('')
Saves previously recorded macro into a script file for future use
-------
Syntax: savemacro [macroname]
- If no macroname is supplied one will be interactively asked
+ If no macroname is supplied one will be interactively asked
'''
saved_ok=0
if self.currentmacro==None:
print 'No macro is being recorded!'
return 0
- if len(macroname)==0:
+ if len(macroname)==0:
macroname=linput.safeinput('Enter new macro name: ')
if len(macroname) == 0:
print 'Invalid name'
macroname.hkm should be present at [hooke]/macros directory
By default the macro will be executed over current curve
- passing 'playlist' word as second argument executes macroname
+ passing 'playlist' word as second argument executes macroname
over all curves
- By default curve(s) will be processed silently, passing 'v'
- as second/third argument will print each command that is
+ By default curve(s) will be processed silently, passing 'v'
+ as second/third argument will print each command that is
executed
Note that macros applied to playlists should end by export
print 'Could not find a macro named '+macropath
return 0
txtfile=open(macropath)
- if cycle ==1:
+ if cycle ==1:
#print self.current_list
for item in self.current_list:
self.current=item
- self.do_plot(0)
+ self.do_plot(0)
for command in txtfile:
if verbose==1:
print 'Executing command '+' '.join(testcmd)
self.onecmd(' '.join(testcmd))
-
-
-
-
-
-
-#!/usr/bin/env python
-
'''
massanalysis.py
class massanalysisCommands:
def _plug_init(self):
- self.mass_variables={}
+ self.mass_variables={}
self.interesting_variables=['curve','firstpeak_distance','lastpeak_distance','Npeaks','median_distance','mean_distance']
self._clean_data()
-
+
def _clean_data(self):
for variable in self.interesting_variables:
self.mass_variables[variable]=[]
'''
calculates X distance of a peak from the contact point
'''
- item.identify(self.drivers)
-
+ item.identify(self.drivers)
+
real_positions=[]
cut_index=self.find_contact_point()
-
+
#we assume the first is the plot with the force curve
plot=item.curve.default_plots()[0]
xret=plot.vectors[1][0]
-
+
start_x=xret[cut_index]
-
- real_positions=[abs((xret[index])-(start_x)) for index in locations]
+
+ real_positions=[abs((xret[index])-(start_x)) for index in locations]
#close all open files
item.curve.close_all()
#needed to avoid *big* memory leaks!
'''
self._clean_data() #if we recall it, clean previous data!
min_deviation=self.convfilt_config['mindeviation']
-
-
+
+
c=0
for item in self.current_list:
try:
peak_location,peak_size=self.exec_has_peaks(item, min_deviation)
real_positions=self.peak_position_from_contact(item, peak_location)
-
+
self.mass_variables['Npeaks'].append(len(peak_location))
-
+
if len(peak_location) > 1:
self.mass_variables['firstpeak_distance'].append(min(real_positions))
self.mass_variables['lastpeak_distance'].append(max(real_positions))
-
+
distancepeaks=[]
for index in range(len(real_positions)-1):
distancepeaks.append(real_positions[index+1]-real_positions[index])
else:
self.mass_variables['firstpeak_distance'].append(0)
- self.mass_variables['lastpeak_distance'].append(0)
-
+ self.mass_variables['lastpeak_distance'].append(0)
+
if len(peak_location) > 2:
self.mass_variables['median_distance'].append(np.median(distancepeaks))
- self.mass_variables['mean_distance'].append(np.mean(distancepeaks))
+ self.mass_variables['mean_distance'].append(np.mean(distancepeaks))
else:
self.mass_variables['median_distance'].append(0)
- self.mass_variables['mean_distance'].append(0)
-
+ self.mass_variables['mean_distance'].append(0)
+
print 'curve',c
except SyntaxError:
print 'curve',c,'not mapped'
pass
-
+
c+=1
def do_plotmap(self,args):
print 'Give me two arguments between those:'
print self.interesting_variables
return
-
+
scattermap=lhc.PlotObject()
scattermap.vectors=[[]]
scattermap.vectors[0].append(x)
scattermap.vectors[0].append(y)
-
+
scattermap.units=[args[0],args[1]]
scattermap.styles=['scatter']
scattermap.destination=1
-
+
self._send_plot([scattermap])
-
+
def do_savemaps(self,args):
'''
args=filename
'''
-
+
'''
def csv_write_cols(data, f):
-
+
#from Bruno Desthuillers on comp.lang.python
-
+
writer = csv.writer(f)
keys = data.keys()
writer.writerow(dict(zip(keys,keys)))
for row in zip(*data.values()):
- writer.writerow(dict(zip(keys, row)))
+ writer.writerow(dict(zip(keys, row)))
'''
-
+
f=open(args,'wb')
lh.csv_write_dictionary(f,self.mass_variables)
f.close()
-
-
\ No newline at end of file
class multidistanceCommands:
-
+
def do_multidistance(self,args):
'''
MULTIDISTANCE
of an existing file, autopeak will resume it and append measurements to it. If you are giving
a new filename, it will create the file and append to it until you close Hooke.
You can also define a minimun deviation of the peaks.
-
+
Syntax:
multidistance [deviation]
deviation = number of times the convolution signal is above the noise absolute deviation.
pk_loc,peak_size=self.has_peaks(defplot, abs_devs)
return pk_loc, peak_size
-
+
noflatten=False
peaks_location, peak_size=find_current_peaks(noflatten, args)
-
+
#if no peaks, we have nothing to plot. exit.
if len(peaks_location)==0:
return
-
+
#otherwise, we plot the peak locations.
xplotted_ret=self.plots[0].vectors[1][0]
yplotted_ret=self.plots[0].vectors[1][1]
xgood=[xplotted_ret[index] for index in peaks_location]
ygood=[yplotted_ret[index] for index in peaks_location]
-
+
recplot=self._get_displayed_plot()
recplot.vectors.append([xgood,ygood])
if recplot.styles==[]:
if exclude_raw=='N':
print 'Discarded.'
return
-
+
if not exclude_raw=='':
exclude=exclude_raw.split(',')
#we convert in numbers the input
peaks_location= peaks_location[0:new_a]+peaks_location[new_a+1:]
peak_size= peak_size[0:new_a]+peak_size[new_a+1:]
count+=1
-
+
#we calculate the distance vector
dist=[]
for i in range(len(peaks_location)-1):
dist.append(xplotted_ret[peaks_location[i]]-xplotted_ret[peaks_location[i+1]])
-
-
+
+
if self.autofile=='':
print 'Not saved.'
return
-
+
if not os.path.exists(self.autofile):
f=open(self.autofile,'w+')
f.write('Analysis started '+time.asctime()+'\n')
f.write(str(o))
f.write("\n")
f.close()
-
+
print 'Saving...'
f=open(self.autofile,'a+')
-
+
f.write(self.current.path+'\n')
for i in dist:
f.write(";")
f.write(str(i))
- f.write("\n")
+ f.write("\n")
f.close()
-#!/usr/bin/env python
-
from mdp import pca
from libhooke import WX_GOOD, ClickedPoint
import wxversion
class pclusterCommands:
-
+
def _plug_init(self):
self.clustplot1=None
self.clustplot2=None
self.my_work_dir = os.getcwd()+slash+pclus_dir+slash
self.my_curr_dir = os.path.basename(os.getcwd())
os.mkdir(self.my_work_dir)
-
+
#--Custom persistent length
pl_value=None
for arg in args.split():
pl_value=float(pl_expression[1]) #actual value
else:
pl_value=None
-
+
#configuration variables
min_npks = self.convfilt_config['minpeaks']
min_deviation = self.convfilt_config['mindeviation']
-
+
pclust_filename = "automeasure_"+self.my_curr_dir+"_blind"+blindw+".txt" #raw_input('Automeasure filename? ')
realclust_filename = "coordinate_"+self.my_curr_dir+"_blind"+blindw+".txt" #raw_input('Coordinates filename? ')
peackforce_filename = "peakforce_"+self.my_curr_dir+"_blind"+blindw+".txt" #raw_input('Peacks and Forces filename? ')
-
+
f=open(self.my_work_dir+pclust_filename,'w+')
f.write('Analysis started '+time.asctime()+'\n')
f.write('----------------------------------------\n')
f.write('; Contour length (nm) ; Persistence length (nm) ; Max.Force (pN) ; Slope (N/m) ; Sigma contour (nm) ; Sigma persistence (nm)\n')
f.close()
-
+
f=open(self.my_work_dir+realclust_filename,'w+')
f.write('Analysis started '+time.asctime()+'\n')
f.write('----------------------------------------\n')
f.write('; Peak number ; Mean delta (nm) ; Median delta (nm) ; Mean force (pN) ; Median force (pN) ; First peak length (nm) ; Last peak length (nm) ; Max force (pN) ; Min force (pN) ; Max delta (nm) ; Min delta (nm) ; Peaks Diff\n')
f.close()
-
+
f=open(self.my_work_dir+peackforce_filename,'w+')
f.write('Analysis started '+time.asctime()+'\n')
f.write('----------------------------------------\n')
f.write('; Peak number ; 1 peak Length (nm) ; 1 peak Force (pN) ; 2 peak Length (nm) ; 2 peak Force (pN) ; 3 peak Length (nm) ; 3 peak Force (pN) ; 4 peak Length (nm) ; 4 peak Force (pN) ; 5 peak Length (nm) ; 5 peak Force (pN) ; 6 peak Length (nm) ; 6 peak Force (pN) ; 7 peak Length (nm) ; 7 peak Force (pN) ; 8 peak Length (nm) ; 8 peak Force (pN)\n')
f.close()
-
+
# ------ FUNCTION ------
def fit_interval_nm(start_index,plot,nm,backwards):
'''
'''
whatset=1 #FIXME: should be decidable
x_vect=plot.vectors[1][0]
-
+
c=0
i=start_index
start=x_vect[start_index]
fit_errors [6.5817195369767644e-010, 2.4415923138871498e-011]
'''
fit_points=int(self.config['auto_fit_points']) # number of points to fit before the peak maximum <50>
-
+
T=self.config['temperature'] #temperature of the system in kelvins. By default it is 293 K. <301.0>
cindex=self.find_contact_point(itplot[0]) #Automatically find contact point <158, libhooke.ClickedPoint>
contact_point=self._clickize(itplot[0].vectors[1][0], itplot[0].vectors[1][1], cindex)
def features_peaks(itplot, peak, fit_points, contact_point, pl_value, T, cindex, avg):
'''
- calculate informations for each peak and add they in
+ calculate informations for each peak and add they in
c_lengths, p_lengths, sigma_c_lengths, sigma_p_lengths, forces, slopes
'''
c_leng=None
sigma_p_leng=None
force=None
slope=None
-
+
delta_force=10
slope_span=int(self.config['auto_slope_span'])
-
+
peak_point=self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],peak)
other_fit_point=self._clickize(itplot[0].vectors[1][0],itplot[0].vectors[1][1],peak-fit_points)
-
+
points=[contact_point, peak_point, other_fit_point]
-
+
params, yfit, xfit, fit_errors = self.wlc_fit(points, itplot[0].vectors[1][0], itplot[0].vectors[1][1], pl_value, T, return_errors=True)
-
+
#Measure forces
delta_to_measure=itplot[0].vectors[1][1][peak-delta_force:peak+delta_force]
y=min(delta_to_measure)
#check if persistent length makes sense. otherwise, discard peak.
if p_leng>self.config['auto_min_p'] and p_leng<self.config['auto_max_p']:
'''
- p_lengths.append(p_leng)
+ p_lengths.append(p_leng)
c_lengths.append(params[0]*(1.0e+9))
sigma_c_lengths.append(fit_errors[0]*(1.0e+9))
sigma_p_lengths.append(fit_errors[1]*(1.0e+9))
forces.append(abs(y-avg)*(1.0e+12))
- slopes.append(slope)
+ slopes.append(slope)
'''
c_leng=params[0]*(1.0e+9)
sigma_c_leng=fit_errors[0]*(1.0e+9)
itplot[0]=flatten(itplot[0], item, customvalue=1)
try:
peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
- except:
+ except:
#We have troubles with exec_has_peaks (bad curve, whatever).
#Print info and go to next cycle.
print 'Cannot process ',item.path
- continue
+ continue
if len(peak_location)==0:
print 'No peaks!'
continue
fit_points, contact_point, pl_value, T, cindex, avg = plot_informations(itplot,pl_value)
-
+
print '\n\nCurve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'
#initialize output data vectors
if len(vect)==0:
for i in range(len(c_lengths)):
vect.append(0)
-
+
print 'Measurements for all peaks detected:'
print 'contour (nm)', c_lengths
print 'sigma contour (nm)',sigma_c_lengths
print 'sigma p (nm)',sigma_p_lengths
print 'forces (pN)',forces
print 'slopes (N/m)',slopes
-
+
'''
write automeasure text file
'''
for i in range(len(c_lengths)):
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')
f.close()
-
+
peak_number=len(c_lengths)
-
+
'''
write peackforce text file
'''
peackforce_info = peackforce_info + ' ; ' + str(c_lengths[i]) + ' ; ' + str(forces[i])
f.write(' ; '+str(peak_number)+peackforce_info+'\n')
f.close()
-
+
'''
calculate clustering coordinates
'''
deltas=[]
for i in range(len(c_lengths)-1):
deltas.append(c_lengths[i+1]-c_lengths[i])
-
+
delta_mean=np.mean(deltas)
delta_median=np.median(deltas)
-
+
force_mean=np.mean(forces)
force_median=np.median(forces)
-
+
first_peak_cl=c_lengths[0]
last_peak_cl=c_lengths[-1]
-
+
max_force=max(forces[:-1])
min_force=min(forces)
-
+
max_delta=max(deltas)
min_delta=min(deltas)
-
+
delta_stdev=np.std(deltas)
forces_stdev=np.std(forces[:-1])
-
+
peaks_diff=(last_peak_cl-first_peak_cl)/peak_number
-
+
print 'Coordinates'
print 'Peaks',peak_number
print 'Mean delta',delta_mean
print 'Delta stdev',delta_stdev
print 'Forces stdev',forces_stdev
print 'Peaks difference',peaks_diff
-
+
'''
write clustering coordinates
'''
' ; '+str(peaks_diff)+ # 12
'\n')
f.close()
-
+
# start PCA
self.do_pca(pclus_dir+"/"+realclust_filename)
-
-
+
+
def do_pca(self,args):
'''
PCA -> "pca gaeta_coor_blind50.txt 1,3,6"
Automatically measures pca from coordinates filename and shows two interactives plots
- With the second argument (arbitrary) you can select the columns and the multiplier factor
+ With the second argument (arbitrary) you can select the columns and the multiplier factor
to use for the pca (for es "1,3*50,6,8x10,9"). Dont use spaces. "*" or "x" are the same thing.
Without second argument it reads pca_config.txt file
(c)Paolo Pancaldi, Massimo Sandal 2009
'''
-
+
# reads the columns of pca
if self.config['hookedir'][0]=='/':
slash='/' #a Unix or Unix-like system
conf=open(self.my_hooke_dir+"pca_config.txt")
config = conf.readlines()
conf.close()
-
+
self.plot_myCoord = [] # tiene le coordinate prese direttamente dal file creato con pCluster
self.plot_origCoord = [] # tiene le coordinate solo delle colonne scelte e moltiplicate per i valori scelti
self.plot_pcaCoord = [] # tiene le due colonne della pca
self.plot_NewPcaCoord = [] # tiene le due colonne della pca filtrate per densita
self.plot_NewPcaCoordTr=[] # tiene le due colonne della pca trasposta filtrate per densita
plot_path_temp = ""
-
- # prende in inpunt un arg (nome del file)
+
+ # prende in inpunt un arg (nome del file)
# e il secondo le colonne su cui lavorare (e' arbitrario, riceve x es "1,2,3")
arg = args.split(" ")
if arg[0]==args:
else:
file_name=arg[0]
config[0] = arg[1]
-
+
# creo l'array "plot_myCoord" con tutte le coordinate dei plots
# e l'array plot_paths con tutti i percorsi dei plots
nPlotTot = -3 #tolgo le prime 3 righe iniziali del file
if row[0]==" " and row.find('nan')==-1 and row.find("-1.#IND")==-1:
row = row[row.index(";",2)+2:].split(" ; ") # non considero la prima colonna col #picchi
row = [float(i) for i in row]
-
+
#0:Mean delta, 1:Median delta, 2:Mean force, 3:Median force, 4:First peak length, 5:Last peak length
#6:Max delta 7:Min delta 8:Max force 9:Min force 10:Std delta 11:Std force
if (row[0]<500 and row[1]<500 and row[2]<500 and row[3]<500 and row[4]<500 and row[5]<500 and row[6]<500 and row[7]<500 and row[8]<500 and row[9]<500 and row[10]<500 and row[11]<500):
self.plot_myCoord.append(row)
self.plot_paths.append(plot_path_temp)
f.close()
-
- # creo l'array con alcune colonne e pure moltiplicate
+
+ # creo l'array con alcune colonne e pure moltiplicate
for row in self.plot_myCoord:
res=[]
for cols in config[0].split(","):
molt = 1
res.append(row[col]*molt)
self.plot_origCoord.append(res)
-
+
# array convert, calculate PCA, transpose
self.plot_origCoord = np.array(self.plot_origCoord,dtype='float')
#print self.plot_origCoord.shape
self.plot_pcaCoordTr = np.transpose(self.plot_pcaCoord)
pca_X=np.array(self.plot_pcaCoordTr[0],dtype='float')
pca_Y=np.array(self.plot_pcaCoordTr[1],dtype='float')
-
+
'''
# Start section of testing with good plots # 4 TESTING!
Xsyn_1=[]
- Ysyn_1=[]
+ Ysyn_1=[]
Xgb1_1=[]
Ygb1_1=[]
Xbad_1=[]
goodnames=goodnamefile.readlines()
nPlotGood = len(goodnames)-2 #tolgo prima e ultima riga
goodnames=[i.split()[0] for i in goodnames[1:]]
-
+
for index in range(len(self.plot_paths)):
if self.plot_paths[index][:-1] in goodnames:
Xsyn_1.append(pca_X[index])
Ybad_1.append(pca_Y[index])
# Stop section of testing with good plots # 4 TESTING!
'''
-
+
# print first plot
clustplot1=lhc.PlotObject()
clustplot1.add_set(pca_X,pca_Y)
clustplot1.destination=0
self._send_plot([clustplot1])
self.clustplot1=clustplot1
-
+
# density and filer estimation
kernel = sp.stats.kde.gaussian_kde(sp.c_[pca_X,pca_Y].T)
tallest = 0
axis_X = np.linspace(xmin,xmax,num=100)
axis_Y = np.linspace(ymin,ymax,num=100)
'''
-
+
# density filtering:
# tramite "kernel.evaluate" trovo lo score (altezza) di ogni coordinata e decido se mantenerla o no
filtered_pca_X = []
filtered_PcaCoordTr.append(filtered_pca_X)
filtered_PcaCoordTr.append(filtered_pca_Y)
filtered_PcaCoord = np.transpose(filtered_PcaCoordTr)
-
+
# creo i due array "plot_FiltOrigCoord" e "plot_FiltPaths" contenenti solo i dati filtrati con alta densita
for index in range(len(self.plot_pcaCoord)):
if self.plot_pcaCoord[index] in filtered_PcaCoord:
self.plot_FiltOrigCoord.append(self.plot_myCoord[index])
self.plot_FiltPaths.append(self.plot_paths[index])
-
+
'''
# START PCA#2: USELESS!!!
-
+
# creo l array con alcune colonne e pure moltiplicate
temp_coord = []
for row in self.plot_FiltOrigCoord:
res.append(row[col]*molt)
temp_coord.append(res)
self.plot_FiltOrigCoord = temp_coord
-
+
# ricalcolo la PCA: array convert, calculate PCA, transpose
self.plot_FiltOrigCoord = np.array(self.plot_FiltOrigCoord,dtype='float')
#print self.plot_FiltOrigCoord.shape
self.plot_NewPcaCoordTr = np.transpose(self.plot_NewPcaCoord)
pca_X2=np.array(self.plot_NewPcaCoordTr[0],dtype='float')
pca_Y2=np.array(self.plot_NewPcaCoordTr[1],dtype='float')
-
+
# Start section of testing with good plots # 4 TESTING!
Xsyn_2=[]
Ysyn_2=[]
else:
Xbad_2.append(pca_X2[index])
Ybad_2.append(pca_Y2[index])
-
+
# print second plot
clustplot2=lhc.PlotObject()
#clustplot2.add_set(pca_X2,pca_Y2)
self._send_plot([clustplot2])
self.clustplot2=clustplot2
'''
-
+
# PRINT density plot
clustplot2=lhc.PlotObject()
clustplot2.add_set(filtered_pca_X,filtered_pca_Y)
clustplot2.destination=1
self._send_plot([clustplot2])
self.clustplot2=clustplot2
-
+
# printing results
config_pca1 = config[0].replace("*", "x").rstrip("\n")
config_pca2 = config[2].replace("*", "x").rstrip("\n")
print "Piu alta: ", tallest
#print "- FILTRO 3: 2'PCA:"+config_pca2
print ""
-
+
# -- exporting coordinates and plot of PCA in debug mode! --
if config[3].find("true")!=-1:
#1' PCA: save plot and build coordinate s file
f.write (" ; " + str(cel))
f.write ("\n")
f.close()
-
+
# pCLUSTER SAVING!!!
import shutil
pcl_name = file_name.replace("coordinate_", "goodplots_").replace('.txt','_'+config_pca1+'_'+str(my_filter).replace(".",","))
f.write (myfile+"\n")
shutil.copy2(myfile, pcl_name)
f.close()
-
-
+
+
def do_multipca(self,args):
'''
MULTIPCA -> "multipca gaeta_coor_blind50.txt 3"
- Automatically multiply the column suggest in second argument for value between 1-100 (step of 2),
+ Automatically multiply the column suggest in second argument for value between 1-100 (step of 2),
measures pca from coordinates filename and save the png plots.
(c)Paolo Pancaldi, Massimo Sandal 2009
'''
for j in range(1, 13):
if i!=j:
self.do_pca(file_name + " " + str(i) + "," + str(j))
-
+
def do_triplepca(self,args):
'''
TRIPLEPCA -> "triplepca gaeta_coor_blind50.txt"
'''
It returns id, coordinates and file name of a clicked dot on a PCA graphic
'''
-
+
self._send_plot([self.clustplot1]) #quick workaround for BAD problems in the GUI
print 'Click point'
point = self._measure_N_points(N=1, whatset=0)
print "coord: " + str(dot_coord)
self.do_genlist(str(plot_file))
#self.do_jump(str(plot_file))
-
+
# indea iniziata e messa da parte...
def do_peakforce(self, args):
'''
Automatically measures peack and force plots
(c)Paolo Pancaldi, Massimo Sandal 2009
'''
-
+
# prende in inpunt un arg (nome del file)
file_name=args
f=open(file_name)
-
+
# scrivo un file temp
g = open('_prove.txt','w')
-
+
plot_file = '';
rows = f.readlines()
for row in rows:
g.write(writeme)
g.close()
f.close()
-
-
\ No newline at end of file
-#!/usr/bin/env python
-
'''
PROCPLOTS
Processed plots plugin for force curves.
import copy
class procplotsCommands:
-
+
def _plug_init(self):
pass
-
+
def do_derivplot(self,args):
'''
DERIVPLOT
Plots the derivate (actually, the discrete differentiation) of the current force curve
---------
Syntax: derivplot
- '''
- dplot=self.derivplot_curves()
- plot_graph=self.list_of_events['plot_graph']
+ '''
+ dplot=self.derivplot_curves()
+ plot_graph=self.list_of_events['plot_graph']
wx.PostEvent(self.frame,plot_graph(plots=[dplot]))
-
+
def derivplot_curves(self):
'''
do_derivplot helper function
-
+
create derivate plot curves for force curves.
- '''
+ '''
dplot=lhc.PlotObject()
- dplot.vectors=[]
-
+ dplot.vectors=[]
+
for vector in self.plots[0].vectors:
dplot.vectors.append([])
dplot.vectors[-1].append(vector[0][:-1])
dplot.vectors[-1].append(np.diff(vector[1]))
-
+
dplot.destination=1
dplot.units=self.plots[0].units
-
- return dplot
-
+
+ return dplot
+
def do_subtplot(self,args):
'''
SUBTPLOT
Syntax: subtplot
'''
#FIXME: sub_filter and sub_order must be args
-
+
if len(self.plots[0].vectors) != 2:
print 'This command only works on a curve with two different plots.'
pass
-
+
outplot=self.subtract_curves(sub_order=1)
-
- plot_graph=self.list_of_events['plot_graph']
+
+ plot_graph=self.list_of_events['plot_graph']
wx.PostEvent(self.frame,plot_graph(plots=[outplot]))
-
+
def subtract_curves(self, sub_order):
'''
subtracts the extension from the retraction
- '''
+ '''
xext=self.plots[0].vectors[0][0]
yext=self.plots[0].vectors[0][1]
xret=self.plots[0].vectors[1][0]
yret=self.plots[0].vectors[1][1]
-
+
#we want the same number of points
maxpoints_tot=min(len(xext),len(xret))
xext=xext[0:maxpoints_tot]
yext=yext[0:maxpoints_tot]
xret=xret[0:maxpoints_tot]
yret=yret[0:maxpoints_tot]
-
+
if sub_order:
ydiff=[yretval-yextval for yretval,yextval in zip(yret,yext)]
else: #reverse subtraction (not sure it's useful, but...)
- ydiff=[yextval-yretval for yextval,yretval in zip(yext,yret)]
-
+ ydiff=[yextval-yretval for yextval,yretval in zip(yext,yret)]
+
outplot=copy.deepcopy(self.plots[0])
outplot.vectors[0][0], outplot.vectors[1][0] = xext,xret #FIXME: if I use xret, it is not correct!
outplot.vectors[1][1]=ydiff
outplot.vectors[0][1]=[0 for item in outplot.vectors[1][0]]
-
+
return outplot
median_filter=customvalue
else:
median_filter=self.config['medfilt']
-
+
if median_filter==0:
return plot
-
+
if float(median_filter)/2 == int(median_filter)/2:
median_filter+=1
-
+
nplots=len(plot.vectors)
c=0
while c<nplots:
plot.vectors[c][1]=scipy.signal.medfilt(plot.vectors[c][1],median_filter)
c+=1
-
+
return plot
-
+
def plotmanip_correct(self, plot, current, customvalue=None):
'''
- the deflection() vector is as long as the X of extension + the X of retraction
- plot.vectors[0][0] is the X of extension curve
- plot.vectors[1][0] is the X of retraction curve
-
+
FIXME: both this method and the picoforce driver have to be updated, deflection() must return
a more senseful data structure!
'''
#use only for force spectroscopy experiments!
if current.curve.experiment != 'smfs':
return plot
-
+
if customvalue != None:
execute_me=customvalue
else:
execute_me=self.config['correct']
if not execute_me:
return plot
-
+
defl_ext,defl_ret=current.curve.deflection()
#halflen=len(deflall)/2
-
+
plot.vectors[0][0]=[(zpoint-deflpoint) for zpoint,deflpoint in zip(plot.vectors[0][0],defl_ext)]
plot.vectors[1][0]=[(zpoint-deflpoint) for zpoint,deflpoint in zip(plot.vectors[1][0],defl_ret)]
return plot
-
+
'''
def plotmanip_detriggerize(self, plot, current, customvalue=None):
#DEPRECATED
if self.config['detrigger']==0:
return plot
-
+
cutindex=2
startvalue=plot.vectors[0][0][0]
-
- for index in range(len(plot.vectors[0][0])-1,2,-2):
+
+ for index in range(len(plot.vectors[0][0])-1,2,-2):
if plot.vectors[0][0][index]>startvalue:
cutindex=index
else:
plot.vectors[0][0]=plot.vectors[0][0][:cutindex]
plot.vectors[0][1]=plot.vectors[0][1][:cutindex]
-
+
return plot
- '''
-
-
-
+ '''
+
+
+
#FFT---------------------------
def fft_plot(self, vector):
'''
'''
fftplot=lhc.PlotObject()
fftplot.vectors=[[]]
-
+
fftlen=len(vector)/2 #need just 1/2 of length
- fftplot.vectors[-1].append(np.arange(1,fftlen).tolist())
-
+ fftplot.vectors[-1].append(np.arange(1,fftlen).tolist())
+
try:
fftplot.vectors[-1].append(abs(np.fft(vector)[1:fftlen]).tolist())
except TypeError: #we take care of newer NumPy (1.0.x)
fftplot.vectors[-1].append(abs(np.fft.fft(vector)[1:fftlen]).tolist())
-
-
+
+
fftplot.destination=1
-
-
+
+
return fftplot
-
-
+
+
def do_fft(self,args):
'''
FFT
Plots the fast Fourier transform of the selected plot
---
Syntax: fft [top,bottom] [select] [0,1...]
-
+
By default, fft performs the Fourier transform on all the 0-th data set on the
top plot.
-
+
[top,bottom]: which plot is the data set to fft (default: top)
[select]: you pick up two points on the plot and fft only the segment between
[0,1,...]: which data set on the selected plot you want to fft (default: 0)
'''
-
+
#args parsing
#whatplot = plot to fft
#whatset = set to fft in the plot
whatset=int(arg)
except ValueError:
pass
-
+
if select:
points=self._measure_N_points(N=2, whatset=whatset)
boundaries=[points[0].index, points[1].index]
y_to_fft=self.plots[whatplot].vectors[whatset][1][boundaries[0]:boundaries[1]] #y
else:
y_to_fft=self.plots[whatplot].vectors[whatset][1] #y
-
- fftplot=self.fft_plot(y_to_fft)
+
+ fftplot=self.fft_plot(y_to_fft)
fftplot.units=['frequency', 'power']
- plot_graph=self.list_of_events['plot_graph']
+ plot_graph=self.list_of_events['plot_graph']
wx.PostEvent(self.frame,plot_graph(plots=[fftplot]))
-
\ No newline at end of file
-#!/usr/bin/env python
-
from libhooke import WX_GOOD
import wxversion
wxversion.select(WX_GOOD)
from numpy import arange, mean
class superimposeCommands:
-
+
def _plug_init(self):
self.imposed=[]
-
+
def do_selimpose(self,args):
'''
SELIMPOSE (superimpose.py plugin)
Hand-selects the curve portion to superimpose
'''
#fixme: set to superimpose should be in args
-
+
if args=='clear':
self.imposed=[]
return
-
+
current_set=1
-
+
points=self._measure_two_points()
boundaries=[points[0].index, points[1].index]
boundaries.sort()
-
+
theplot=self.plots[0]
#append the selected section
self.imposed.append([])
self.imposed[-1].append(theplot.vectors[1][0][boundaries[0]:boundaries[1]]) #x
self.imposed[-1].append(theplot.vectors[1][1][boundaries[0]:boundaries[1]]) #y
-
+
#align X first point
self.imposed[-1][0] = [item-self.imposed[-1][0][0] for item in self.imposed[-1][0]]
#align Y first point
self.imposed[-1][1] = [item-self.imposed[-1][1][0] for item in self.imposed[-1][1]]
-
+
def do_plotimpose(self,args):
'''
PLOTIMPOSE (sumperimpose.py plugin)
imposed_object=lhc.PlotObject()
imposed_object.vectors=self.imposed
print 'Plotting',len(imposed_object.vectors),'imposed curves'
-
+
imposed_object.normalize_vectors()
-
+
imposed_object.units=self.plots[0].units
imposed_object.title='Imposed curves'
imposed_object.destination=1
-
- plot_graph=self.list_of_events['plot_graph']
+
+ plot_graph=self.list_of_events['plot_graph']
PostEvent(self.frame,plot_graph(plots=[imposed_object]))
-
+
def do_plotavgimpose(self,args):
'''
PLOTAVGIMPOSE (superimpose.py plugin)
min_x=[]
for curve in self.imposed:
min_x.append(min(curve[0]))
-
+
#find minimum extension
min_ext_limit=max(min_x)
-
+
x_avg=arange(step,min_ext_limit,step)
y_avg=[]
for value in x_avg:
for curve in self.imposed:
for xvalue, yvalue in zip(curve[0],curve[1]):
if xvalue >= (value+step) and xvalue <= (value-step):
- to_avg.append(yvalue)
+ to_avg.append(yvalue)
y_avg.append(mean(to_avg))
-
+
print 'len x len y'
print len(x_avg), len(y_avg)
print y_avg
-
+
avg_object=lhc.PlotObject()
avg_object.vectors=[[x_avg, y_avg]]
avg_object.normalize_vectors()
avg_object.units=self.plots[0].units
avg_object.title="Average curve"
avg_object.destination=1
-
- plot_graph=self.list_of_events['plot_graph']
- PostEvent(self.frame,plot_graph(plots=[avg_object]))
\ No newline at end of file
+
+ plot_graph=self.list_of_events['plot_graph']
+ PostEvent(self.frame,plot_graph(plots=[avg_object]))
-#!/usr/bin/env python
-
'''
TUTORIAL PLUGIN FOR HOOKE
'''
Here we define the class containing all the Hooke commands we want to define
in the plugin.
-
+
Notice the class name!!
The syntax is filenameCommands. That is, if your plugin is pluggy.py, your class
name is pluggyCommands.
-
+
Otherwise, the class will be ignored by Hooke.
- '''
-
+ '''
+
def _plug_init(self):
'''
This is the plugin initialization.
If there is something you need to do when Hooke starts, code it in this function.
'''
print 'I am the Tutorial plugin initialization!'
-
+
#Here we initialize a local configuration variable; see plotmanip_absvalue() for explanation.
- self.config['tutorial_absvalue']=0
+ self.config['tutorial_absvalue']=0
pass
-
+
def do_nothing(self,args):
'''
This is a boring but working example of an actual Hooke command.
A Hooke command is a function of the xxxxCommands class, which is ALWAYS defined
this way:
-
+
def do_nameofcommand(self,args)
-
+
*do_ is needed to make Hooke understand this function is a command
*nameofcommand is how the command will be called in the Hooke command line.
*self is, well, self
*args is ALWAYS needed (otherwise Hooke will crash executing the command). We will see
later what args is.
-
+
Note that if you now start Hooke with this plugin activated and you type in the Hooke command
line "help nothing" you will see this very text as output. So the help of a command is a
string comment below the function definition, like this one.
-
+
Commands usually return None.
'''
print 'I am a Hooke command. I do nothing.'
-
+
def do_printargs(self,args):
'''
This command prints the args you give to it.
args is always a string, that contains everything you write after the command.
So if you issue "mycommand blah blah 12345" args is "blah blah 12345".
-
+
Again, args is needed in the definition even if your command does not use it.
'''
print 'You gave me those args: '+args
-
+
def help_tutorial(self):
'''
- This is a help function.
+ This is a help function.
If you want a help function for something that is not a command, you can write a help
function like this. Calling "help tutorial" will execute this function.
'''
print 'You called help_tutorial()'
-
+
def do_environment(self,args):
'''
This plugin contains a panoramic of the Hooke command line environment variables,
and prints their current value.
'''
-
+
'''self.current_list
TYPE: [ libhookecurve.HookeCurve ], len=variable
contains the actual playlist of Hooke curve objects.
'''
print 'current_list length:',len(self.current_list)
print 'current_list 0th:',self.current_list[0]
-
+
'''self.pointer
TYPE: int
contains the index of
the current curve in the playlist
'''
print 'pointer: ',self.pointer
-
+
'''self.current
TYPE: libhookecurve.HookeCurve
contains the current curve displayed.
We will see later how it works.
'''
print 'current:',self.current
-
+
'''self.plots
TYPE: [ libhookecurve.PlotObject ], len=1,2
contains the current default plots.
- Each PlotObject contains all info needed to display
+ Each PlotObject contains all info needed to display
the plot: apart from the data vectors, the title, destination
etc.
Usually self.plots[0] is the default topmost plot, self.plots[1] is the
accessory bottom plot.
'''
print 'plots:',self.plots
-
+
'''self.config
TYPE: { string:anything }
contains the current Hooke configuration variables, in form of a dictionary.
'''
print 'config:',self.config
-
+
'''self.plotmanip
TYPE: [ function ]
Contains the ordered plot manipulation functions.
*YOU SHOULD NEVER MODIFY THAT*
'''
print 'plotmanip: ',self.plotmanip
-
+
'''self.drivers
TYPE: [ class ]
Contains the plot reading drivers.
*YOU SHOULD NEVER MODIFY THAT*
'''
print 'drivers: ',self.drivers
-
+
'''self.frame
TYPE: wx.Frame
Contains the wx Frame of the GUI.
***NEVER, EVER TOUCH THAT.***
'''
print 'frame: ',self.frame
-
+
'''self.list_of_events
TYPE: { string:wx.Event }
Contains the wx.Events to communicate with the GUI.
to create a GUI plugin.
'''
print 'list of events:',self.list_of_events
-
+
'''self.events_from_gui
TYPE: Queue.Queue
Contains the Queue where data from the GUI is put.
to create a GUI plugin.
'''
print 'events from gui:',self.events_from_gui
-
+
'''self.playlist_saved
TYPE: Int (0/1) ; Boolean
Flag that tells if the playlist has been saved or not.
'''
print 'playlist saved:',self.playlist_saved
-
+
'''self.playlist_name
TYPE: string
Name of current playlist
'''
print 'playlist name:',self.playlist_name
-
+
'''self.notes_saved
TYPE: Int (0/1) ; Boolean
Flag that tells if the playlist has been saved or not.
'''
print 'notes saved:',self.notes_saved
-
+
def do_myfirstplot(self,args):
'''
In this function, we see how to create a PlotObject and send it to the screen.
***Read the code of PlotObject in libhookecurve.py before!***.
'''
-
+
#We generate some interesting data to plot for this example.
xdata1=np.arange(-5,5,0.1)
xdata2=np.arange(-5,5,0.1)
ydata1=[item**2 for item in xdata1]
ydata2=[item**3 for item in xdata2]
-
+
#Create the object.
#The PlotObject class lives in the libhookecurve library.
myplot=lhc.PlotObject()
'''
- The *data* of the plot live in the .vectors list.
-
+ The *data* of the plot live in the .vectors list.
+
plot.vectors is a multidimensional array:
plot.vectors[0]=set1
plot.vectors[1]=set2
plot.vectors[2]=sett3
etc.
-
+
2 curves in a x,y plot are:
[[[x1],[y1]],[[x2],[y2]]]
for example:
x2 = self.vectors[1][0]
y2 = self.vectors[1][1]
'''
- #Pour 0-th dataset into myplot:
+ #Pour 0-th dataset into myplot:
myplot.add_set(xdata1,ydata1)
-
- #Pour 1-st dataset into myplot:
+
+ #Pour 1-st dataset into myplot:
myplot.add_set(xdata2,ydata2)
-
+
#Add units to x and y axes
#units=[string, string]
myplot.units=['x axis','y axis']
-
+
#Where do we want the plot? 0=top, 1=bottom
myplot.destination=1
-
+
'''Send it to the GUI.
Note that you *have* to encapsulate it into a list, so you
have to send [myplot], not simply myplot.
-
+
You can also send more two plots at once
self.send_plot([plot1,plot2])
'''
self._send_plot([myplot])
-
+
def do_myfirstscatter(self,args):
'''
xdata2=np.arange(-5,5,1)
ydata1=[item**2 for item in xdata1]
ydata2=[item**3 for item in xdata2]
-
+
myplot=lhc.PlotObject()
myplot.add_set(xdata1,ydata1)
myplot.add_set(xdata2,ydata2)
-
-
+
+
#Add units to x and y axes
myplot.units=['x axis','y axis']
-
+
#Where do we want the plot? 0=top, 1=bottom
myplot.destination=1
-
+
'''None=standard line plot
'scatter'=scatter plot
By default, the styles attribute is an empty list. If you
want to define a scatter plot, you must define all other
plots as None or 'scatter', depending on what you want.
-
+
Here we define the second set to be plotted as scatter,
and the first to be plotted as line.
-
+
Here we define also the colors to be the default Matplotlib colors
'''
myplot.styles=[None,'scatter']
myplot.colors=[None,None]
self._send_plot([myplot])
-
+
def do_clickaround(self,args):
'''
Here we click two points on the curve and take some parameters from the points
we have clicked.
'''
-
+
'''
points = self._measure_N_points(N=Int, whatset=Int)
*N = number of points to measure(1...n)
'''
points=self._measure_N_points(N=2,whatset=1)
print 'You clicked the following points.'
-
+
'''
These are the absolute coordinates of the
- point clicked.
+ point clicked.
[float, float] = x,y
'''
print 'Absolute coordinates:'
print points[0].absolute_coords
print points[1].absolute_coords
print
-
+
'''
These are the coordinates of the points
clicked, remapped on the graph.
print points[0].graph_coords
print points[1].graph_coords
print
-
+
'''
These are the indexes of the clicked points
on the dataset vector.
print 'Index of points on the graph:'
print points[0].index
print points[1].index
-
-
+
+
def help_thedifferentplots(self):
'''
The *three* different default plots you should be familiar with
in Hooke.
-
+
Each plot contains of course the respective data in their
vectors attribute, so here you learn also which data access for
each situation.
'''
print '''
1. THE RAW, CURRENT PLOTS
-
+
self.current
---
Contains the current libhookecurve.HookeCurve container object.
A HookeCurve object defines only two default attributes:
-
+
* self.current.path = string
The path of the current displayed curve
-
+
* self.current.curve = libhookecurve.Driver
The curve object. This is not only generated by the driver,
this IS a driver instance in itself.
This means that in self.current.curve you can access the
specific driver APIs, if you know them.
-
+
And defines only one method:
* self.current.identify()
Fills in the self.current.curve object.
See in the cycling tutorial.
-
+
*****
The REAL curve data actually lives in:
---
Contains the raw PlotObject-s, as "spitted out" by the driver, without any
intervention.
This is as close to the raw data as Hooke gets.
-
+
One or two plots can be spit out; they are always enclosed in a list.
*****
-
+
Methods of self.current.curve are:
---
-
+
* self.current.curve.is_me()
(Used by identify() only.)
-
+
* self.current.curve.close_all()
Closes all driver open files; see the cycling tutorial.
'''
-
+
print '''
2. THE PROCESSED, DEFAULT PLOT
-
+
The plot that is spitted out by the driver is *not* the usual default plot
that is displayed by calling "plot" at the Hooke prompt.
-
+
This is because the raw, driver-generated plot is usually *processed* by so called
*plot processing* functions. We will see in the tutorial how to define
- them.
-
+ them.
+
For example, in force spectroscopy force curves, raw data are automatically corrected
for deflection. Other data can be, say, filtered by default.
-
- The default plots are accessible in
+
+ The default plots are accessible in
self.plots = [ libhooke.PlotObject ]
-
+
self.plots[0] is usually the topmost plot
self.plots[1] is usually the bottom plot (if present)
'''
-
+
print '''
3. THE PLOT DISPLAYED RIGHT NOW.
-
+
Sometimes the plots you are displaying *right now* is different from the previous
two. You may have a fit trace, you may have issued some command that spits out
- a custom plot and you want to rework that, whatever.
-
+ a custom plot and you want to rework that, whatever.
+
You can obtain in any moment the plot currently displayed by Hooke by issuing
-
+
PlotObject = self._get_displayed_plot(dest)
* dest = Int (0/1)
dest=0 : top plot
dest=1 : bottom plot
'''
-
-
+
+
def do_cycling(self,args):
'''
Here we cycle through our playlist and print some info on the curves we find.
Cycling through the playlist needs a bit of care to avoid memory leaks and dangling
open files...
-
+
Look at the source code for more information.
'''
-
+
def things_when_cycling(item):
'''
We encapsulate here everything has to open the actual curve file.
By doing it all here, we avoid to do acrobacies when deleting objects etc.
in the main loop: we do the dirty stuff here.
'''
-
+
'''
identify()
-
+
This method looks for the correct driver in self.drivers to use;
and puts the curve content in the .curve attribute.
Basically, until identify() is called, the HookeCurve object
the Hooke plot routine), the HookeCurve object is "filled" with
the actual curve.
'''
-
+
item.identify(self.drivers)
-
+
'''
After the identify(), item.curve contains the curve, and item.curve.default_plots() behaves exactly like
self.current.curve.default_plots() -but for the given item.
'''
itplot=item.curve.default_plots()
-
+
print 'length of X1 vector:',len(itplot[0].vectors[0][0]) #just to show something
-
+
'''
The following three lines are a magic spell you HAVE to do
before closing the function.
item.curve.close_all() #Avoid open files dangling
del item.curve #Avoid memory leaks
del item #Just be paranoid. Don't ask.
-
+
return
-
-
+
+
c=0
for item in self.current_list:
print 'Looking at curve ',c,'of',len(self.current_list)
things_when_cycling(item)
c+=1
-
+
return
-
-
-
+
+
+
def plotmanip_absvalue(self, plot, current, customvalue=None):
'''
This function defines a PLOT MANIPULATOR.
A plot manipulator is a function that takes a plot in input, does something to the plot
and returns the modified plot in output.
The function, once plugged, gets automatically called everytime self.plots is updated
-
+
For example, in force spectroscopy force curves, raw data are automatically corrected
for deflection. Other data can be, say, filtered by default.
-
+
To create and activate a plot manipulator you have to:
* Write a function (like this) which name starts with "plotmanip_" (just like commands
start with "do_")
* The function must return a plot object.
* Add an entry in hooke.conf: if your function is "plotmanip_something" you will have
to add <something/> in the plotmanips section: example
-
+
<plotmanips>
<detriggerize/>
<correct/>
<median/>
- <something/>
+ <something/>
</plotmanips>
-
+
Important: Plot manipulators are *in pipe*: each plot manipulator output becomes the input of the next one.
The order in hooke.conf *is the order* in which plot manipulators are connected, so in the example above
we have:
self.current.curve.default_plots() --> detriggerize --> correct --> median --> something --> self.plots
'''
-
+
'''
Here we see what is in a configuration variable to enable/disable the plot manipulator as user will using
the Hooke "set" command.
'''
if not self.config['tutorial_absvalue']:
return plot
-
+
#We do something to the plot, for demonstration's sake
#If we needed variables, we would have used customvalue.
plot.vectors[0][1]=[abs(i) for i in plot.vectors[0][1]]
plot.vectors[1][1]=[abs(i) for i in plot.vectors[1][1]]
-
+
#Return the plot object.
return plot
-
-
+
+
#TODO IN TUTORIAL:
#how to add lines to an existing plot!!
#peaks
-#!/usr/bin/env python
-
'''
Viewer test case
def do_vwnew(self,args):
#creates a new viewer
self.viewerlist.append(lview.Ascii(self.outlet))
- dt=linput.safeinput('What type of data will this viewer handle? (force/distance/all)',['force', 'distance', 'all'])
+ dt=linput.safeinput('What type of data will this viewer handle? (force/distance/all)',['force', 'distance', 'all'])
#TODO update types, make a list somewhere?
print dt
self.viewerlist[-1].setdtype(dt)
if len(args)==0:
args=0
- self.viewerlist[int(args)].action()
\ No newline at end of file
+ self.viewerlist[int(args)].action()
projects / hooke.git / commitdiff