commit 41b252a3c9e1cab23307b7b4610c7f4e79f5d062
Author: Alexander Minges
Date: Fri Jan 30 03:34:25 2015 +0100
First commit
diff --git a/.eric5project/PyDSF.e4q b/.eric5project/PyDSF.e4q
new file mode 100644
index 0000000..114eba5
--- /dev/null
+++ b/.eric5project/PyDSF.e4q
@@ -0,0 +1,6 @@
+
+
+
+
+
+
diff --git a/.eric5project/PyDSF.e4t b/.eric5project/PyDSF.e4t
new file mode 100644
index 0000000..657aa27
--- /dev/null
+++ b/.eric5project/PyDSF.e4t
@@ -0,0 +1,34 @@
+
+
+
+
+
+
+
+ TODO: not implemented yet
+
+ 2012-10-25, 12:56:20
+
+ ui/main.py
+ 32
+
+
+
+ TODO: the window parameter could be the window itself if a list instead of
+
+ 2012-10-30, 16:01:08
+
+ peakdetect.py
+ 594
+
+
+
+ TODO: not implemented yet
+
+ 2012-10-30, 16:26:45
+
+ ui/mainwindow.py
+ 140
+
+
+
diff --git a/main.py b/main.py
new file mode 100644
index 0000000..4690887
--- /dev/null
+++ b/main.py
@@ -0,0 +1,12 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+
+from PyQt4 import QtCore, QtGui
+from ui.mainwindow import MainWindow
+
+if __name__ == "__main__":
+ import sys
+ app = QtGui.QApplication(sys.argv)
+ ui = MainWindow()
+ ui.show()
+ sys.exit(app.exec_())
diff --git a/peakdetect.py b/peakdetect.py
new file mode 100644
index 0000000..5add7e8
--- /dev/null
+++ b/peakdetect.py
@@ -0,0 +1,735 @@
+#!/usr/bin/env python2
+import numpy as np
+from math import pi, log
+import pylab
+from scipy import fft, ifft
+from scipy.optimize import curve_fit
+
+i = 10000
+x = np.linspace(0, 3.5 * pi, i)
+y = (0.3*np.sin(x) + np.sin(1.3 * x) + 0.9 * np.sin(4.2 * x) + 0.06 *
+ np.random.randn(i))
+
+
+def _datacheck_peakdetect(x_axis, y_axis):
+ if x_axis is None:
+ x_axis = range(len(y_axis))
+
+ if len(y_axis) != len(x_axis):
+ raise (ValueError,
+ 'Input vectors y_axis and x_axis must have same length')
+
+ #needs to be a numpy array
+ y_axis = np.array(y_axis)
+ x_axis = np.array(x_axis)
+ return x_axis, y_axis
+
+def _peakdetect_parabole_fitter(raw_peaks, x_axis, y_axis, points):
+ """
+ Performs the actual parabole fitting for the peakdetect_parabole function.
+
+ keyword arguments:
+ raw_peaks -- A list of either the maximium or the minimum peaks, as given
+ by the peakdetect_zero_crossing function, with index used as x-axis
+ x_axis -- A numpy list of all the x values
+ y_axis -- A numpy list of all the y values
+ points -- How many points around the peak should be used during curve
+ fitting, must be odd.
+
+ return -- A list giving all the peaks and the fitted waveform, format:
+ [[x, y, [fitted_x, fitted_y]]]
+
+ """
+ func = lambda x, k, tau, m: k * ((x - tau) ** 2) + m
+ fitted_peaks = []
+ for peak in raw_peaks:
+ index = peak[0]
+ x_data = x_axis[index - points // 2: index + points // 2 + 1]
+ y_data = y_axis[index - points // 2: index + points // 2 + 1]
+ # get a first approximation of tau (peak position in time)
+ tau = x_axis[index]
+ # get a first approximation of peak amplitude
+ m = peak[1]
+
+ # build list of approximations
+ # k = -m as first approximation?
+ p0 = (-m, tau, m)
+ popt, pcov = curve_fit(func, x_data, y_data, p0)
+ # retrieve tau and m i.e x and y value of peak
+ x, y = popt[1:3]
+
+ # create a high resolution data set for the fitted waveform
+ x2 = np.linspace(x_data[0], x_data[-1], points * 10)
+ y2 = func(x2, *popt)
+
+ fitted_peaks.append([x, y, [x2, y2]])
+
+ return fitted_peaks
+
+
+def peakdetect(y_axis, x_axis = None, lookahead = 300, delta=0):
+ """
+ Converted from/based on a MATLAB script at:
+ http://billauer.co.il/peakdet.html
+
+ function for detecting local maximas and minmias in a signal.
+ Discovers peaks by searching for values which are surrounded by lower
+ or larger values for maximas and minimas respectively
+
+ keyword arguments:
+ y_axis -- A list containg the signal over which to find peaks
+ x_axis -- (optional) A x-axis whose values correspond to the y_axis list
+ and is used in the return to specify the postion of the peaks. If
+ omitted an index of the y_axis is used. (default: None)
+ lookahead -- (optional) distance to look ahead from a peak candidate to
+ determine if it is the actual peak (default: 200)
+ '(sample / period) / f' where '4 >= f >= 1.25' might be a good value
+ delta -- (optional) this specifies a minimum difference between a peak and
+ the following points, before a peak may be considered a peak. Useful
+ to hinder the function from picking up false peaks towards to end of
+ the signal. To work well delta should be set to delta >= RMSnoise * 5.
+ (default: 0)
+ delta function causes a 20% decrease in speed, when omitted
+ Correctly used it can double the speed of the function
+
+ return -- two lists [max_peaks, min_peaks] containing the positive and
+ negative peaks respectively. Each cell of the lists contains a tupple
+ of: (position, peak_value)
+ to get the average peak value do: np.mean(max_peaks, 0)[1] on the
+ results to unpack one of the lists into x, y coordinates do:
+ x, y = zip(*tab)
+ """
+ max_peaks = []
+ min_peaks = []
+ dump = [] #Used to pop the first hit which almost always is false
+
+ # check input data
+ x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)
+ # store data length for later use
+ length = len(y_axis)
+
+
+ #perform some checks
+ if lookahead < 1:
+ raise ValueError, "Lookahead must be '1' or above in value"
+ if not (np.isscalar(delta) and delta >= 0):
+ raise ValueError, "delta must be a positive number"
+
+ #maxima and minima candidates are temporarily stored in
+ #mx and mn respectively
+ mn, mx = np.Inf, -np.Inf
+
+ #Only detect peak if there is 'lookahead' amount of points after it
+ for index, (x, y) in enumerate(zip(x_axis[:-lookahead],
+ y_axis[:-lookahead])):
+ if y > mx:
+ mx = y
+ mxpos = x
+ if y < mn:
+ mn = y
+ mnpos = x
+
+ ####look for max####
+ if y < mx-delta and mx != np.Inf:
+ #Maxima peak candidate found
+ #look ahead in signal to ensure that this is a peak and not jitter
+ if y_axis[index:index+lookahead].max() < mx:
+ max_peaks.append([mxpos, mx])
+ dump.append(True)
+ #set algorithm to only find minima now
+ mx = np.Inf
+ mn = np.Inf
+ if index+lookahead >= length:
+ #end is within lookahead no more peaks can be found
+ break
+ continue
+ #else: #slows shit down this does
+ # mx = ahead
+ # mxpos = x_axis[np.where(y_axis[index:index+lookahead]==mx)]
+
+ ####look for min####
+ if y > mn+delta and mn != -np.Inf:
+ #Minima peak candidate found
+ #look ahead in signal to ensure that this is a peak and not jitter
+ if y_axis[index:index+lookahead].min() > mn:
+ min_peaks.append([mnpos, mn])
+ dump.append(False)
+ #set algorithm to only find maxima now
+ mn = -np.Inf
+ mx = -np.Inf
+ if index+lookahead >= length:
+ #end is within lookahead no more peaks can be found
+ break
+ #else: #slows shit down this does
+ # mn = ahead
+ # mnpos = x_axis[np.where(y_axis[index:index+lookahead]==mn)]
+
+
+ #Remove the false hit on the first value of the y_axis
+ try:
+ if dump[0]:
+ max_peaks.pop(0)
+ else:
+ min_peaks.pop(0)
+ del dump
+ except IndexError:
+ #no peaks were found, should the function return empty lists?
+ pass
+
+ return [max_peaks, min_peaks]
+
+
+def peakdetect_fft(y_axis, x_axis, pad_len = 5):
+ """
+ Performs a FFT calculation on the data and zero-pads the results to
+ increase the time domain resolution after performing the inverse fft and
+ send the data to the 'peakdetect' function for peak
+ detection.
+
+ Omitting the x_axis is forbidden as it would make the resulting x_axis
+ value silly if it was returned as the index 50.234 or similar.
+
+ Will find at least 1 less peak then the 'peakdetect_zero_crossing'
+ function, but should result in a more precise value of the peak as
+ resolution has been increased. Some peaks are lost in an attempt to
+ minimize spectral leakage by calculating the fft between two zero
+ crossings for n amount of signal periods.
+
+ The biggest time eater in this function is the ifft and thereafter it's
+ the 'peakdetect' function which takes only half the time of the ifft.
+ Speed improvementd could include to check if 2**n points could be used for
+ fft and ifft or change the 'peakdetect' to the 'peakdetect_zero_crossing',
+ which is maybe 10 times faster than 'peakdetct'. The pro of 'peakdetect'
+ is that it resutls in one less lost peak. It should also be noted that the
+ time used by the ifft function can change greatly depending on the input.
+
+ keyword arguments:
+ y_axis -- A list containg the signal over which to find peaks
+ x_axis -- A x-axis whose values correspond to the y_axis list and is used
+ in the return to specify the postion of the peaks.
+ pad_len -- (optional) By how many times the time resolution should be
+ increased by, e.g. 1 doubles the resolution. The amount is rounded up
+ to the nearest 2 ** n amount (default: 5)
+
+ return -- two lists [max_peaks, min_peaks] containing the positive and
+ negative peaks respectively. Each cell of the lists contains a tupple
+ of: (position, peak_value)
+ to get the average peak value do: np.mean(max_peaks, 0)[1] on the
+ results to unpack one of the lists into x, y coordinates do:
+ x, y = zip(*tab)
+ """
+ # check input data
+ x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)
+ zero_indices = zero_crossings(y_axis, window = 11)
+ #select a n amount of periods
+ last_indice = - 1 - (1 - len(zero_indices) & 1)
+ # Calculate the fft between the first and last zero crossing
+ # this method could be ignored if the begining and the end of the signal
+ # are discardable as any errors induced from not using whole periods
+ # should mainly manifest in the beginning and the end of the signal, but
+ # not in the rest of the signal
+ fft_data = fft(y_axis[zero_indices[0]:zero_indices[last_indice]])
+ padd = lambda x, c: x[:len(x) // 2] + [0] * c + x[len(x) // 2:]
+ n = lambda x: int(log(x)/log(2)) + 1
+ # padds to 2**n amount of samples
+ fft_padded = padd(list(fft_data), 2 **
+ n(len(fft_data) * pad_len) - len(fft_data))
+
+ # There is amplitude decrease directly proportional to the sample increase
+ sf = len(fft_padded) / float(len(fft_data))
+ # There might be a leakage giving the result an imaginary component
+ # Return only the real component
+ y_axis_ifft = ifft(fft_padded).real * sf #(pad_len + 1)
+ x_axis_ifft = np.linspace(
+ x_axis[zero_indices[0]], x_axis[zero_indices[last_indice]],
+ len(y_axis_ifft))
+ # get the peaks to the interpolated waveform
+ max_peaks, min_peaks = peakdetect(y_axis_ifft, x_axis_ifft, 500,
+ delta = abs(np.diff(y_axis).max() * 2))
+ #max_peaks, min_peaks = peakdetect_zero_crossing(y_axis_ifft, x_axis_ifft)
+
+ # store one 20th of a period as waveform data
+ data_len = int(np.diff(zero_indices).mean()) / 10
+ data_len += 1 - data_len & 1
+
+
+ fitted_wave = []
+ for peaks in [max_peaks, min_peaks]:
+ peak_fit_tmp = []
+ index = 0
+ for peak in peaks:
+ index = np.where(x_axis_ifft[index:]==peak[0])[0][0] + index
+ x_fit_lim = x_axis_ifft[index - data_len // 2:
+ index + data_len // 2 + 1]
+ y_fit_lim = y_axis_ifft[index - data_len // 2:
+ index + data_len // 2 + 1]
+
+ peak_fit_tmp.append([x_fit_lim, y_fit_lim])
+ fitted_wave.append(peak_fit_tmp)
+
+ #pylab.plot(range(len(fft_data)), fft_data)
+ #pylab.show()
+
+ #pylab.plot(x_axis, y_axis)
+ #pylab.hold(True)
+ #pylab.plot(x_axis_ifft, y_axis_ifft)
+ #for max_p in max_peaks:
+ # pylab.plot(max_p[0], max_p[1], 'xr')
+ #pylab.show()
+ return [max_peaks, min_peaks]
+
+
+def peakdetect_parabole(y_axis, x_axis, points = 9):
+ """
+ Function for detecting local maximas and minmias in a signal.
+ Discovers peaks by fitting the model function: y = k (x - tau) ** 2 + m
+ to the peaks. The amount of points used in the fitting is set by the
+ points argument.
+
+ Omitting the x_axis is forbidden as it would make the resulting x_axis
+ value silly if it was returned as index 50.234 or similar.
+
+ will find the same amount of peaks as the 'peakdetect_zero_crossing'
+ function, but might result in a more precise value of the peak.
+
+ keyword arguments:
+ y_axis -- A list containg the signal over which to find peaks
+ x_axis -- A x-axis whose values correspond to the y_axis list and is used
+ in the return to specify the postion of the peaks.
+ points -- (optional) How many points around the peak should be used during
+ curve fitting, must be odd (default: 9)
+
+ return -- two lists [max_peaks, min_peaks] containing the positive and
+ negative peaks respectively. Each cell of the lists contains a list
+ of: (position, peak_value)
+ to get the average peak value do: np.mean(max_peaks, 0)[1] on the
+ results to unpack one of the lists into x, y coordinates do:
+ x, y = zip(*max_peaks)
+ """
+ # check input data
+ x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)
+ # make the points argument odd
+ points += 1 - points % 2
+ #points += 1 - int(points) & 1 slower when int conversion needed
+
+ # get raw peaks
+ max_raw, min_raw = peakdetect_zero_crossing(y_axis)
+
+ # define output variable
+ max_peaks = []
+ min_peaks = []
+
+ max_ = _peakdetect_parabole_fitter(max_raw, x_axis, y_axis, points)
+ min_ = _peakdetect_parabole_fitter(min_raw, x_axis, y_axis, points)
+
+ max_peaks = map(lambda x: [x[0], x[1]], max_)
+ max_fitted = map(lambda x: x[-1], max_)
+ min_peaks = map(lambda x: [x[0], x[1]], min_)
+ min_fitted = map(lambda x: x[-1], min_)
+
+
+ #pylab.plot(x_axis, y_axis)
+ #pylab.hold(True)
+ #for max_p, max_f in zip(max_peaks, max_fitted):
+ # pylab.plot(max_p[0], max_p[1], 'x')
+ # pylab.plot(max_f[0], max_f[1], 'o', markersize = 2)
+ #for min_p, min_f in zip(min_peaks, min_fitted):
+ # pylab.plot(min_p[0], min_p[1], 'x')
+ # pylab.plot(min_f[0], min_f[1], 'o', markersize = 2)
+ #pylab.show()
+
+ return [max_peaks, min_peaks]
+
+
+def peakdetect_sine(y_axis, x_axis, points = 9, lock_frequency = False):
+ """
+ Function for detecting local maximas and minmias in a signal.
+ Discovers peaks by fitting the model function:
+ y = A * sin(2 * pi * f * x - tau) to the peaks. The amount of points used
+ in the fitting is set by the points argument.
+
+ Omitting the x_axis is forbidden as it would make the resulting x_axis
+ value silly if it was returned as index 50.234 or similar.
+
+ will find the same amount of peaks as the 'peakdetect_zero_crossing'
+ function, but might result in a more precise value of the peak.
+
+ The function might have some problems if the sine wave has a
+ non-negligible total angle i.e. a k*x component, as this messes with the
+ internal offset calculation of the peaks, might be fixed by fitting a
+ k * x + m function to the peaks for offset calculation.
+
+ keyword arguments:
+ y_axis -- A list containg the signal over which to find peaks
+ x_axis -- A x-axis whose values correspond to the y_axis list and is used
+ in the return to specify the postion of the peaks.
+ points -- (optional) How many points around the peak should be used during
+ curve fitting, must be odd (default: 9)
+ lock_frequency -- (optional) Specifies if the frequency argument of the
+ model function should be locked to the value calculated from the raw
+ peaks or if optimization process may tinker with it. (default: False)
+
+ return -- two lists [max_peaks, min_peaks] containing the positive and
+ negative peaks respectively. Each cell of the lists contains a tupple
+ of: (position, peak_value)
+ to get the average peak value do: np.mean(max_peaks, 0)[1] on the
+ results to unpack one of the lists into x, y coordinates do:
+ x, y = zip(*tab)
+ """
+ # check input data
+ x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)
+ # make the points argument odd
+ points += 1 - points % 2
+ #points += 1 - int(points) & 1 slower when int conversion needed
+
+ # get raw peaks
+ max_raw, min_raw = peakdetect_zero_crossing(y_axis)
+
+ # define output variable
+ max_peaks = []
+ min_peaks = []
+
+ # get global offset
+ offset = np.mean([np.mean(max_raw, 0)[1], np.mean(min_raw, 0)[1]])
+ # fitting a k * x + m function to the peaks might be better
+ #offset_func = lambda x, k, m: k * x + m
+
+ # calculate an approximate frequenzy of the signal
+ Hz = []
+ for raw in [max_raw, min_raw]:
+ if len(raw) > 1:
+ peak_pos = [x_axis[index] for index in zip(*raw)[0]]
+ Hz.append(np.mean(np.diff(peak_pos)))
+ Hz = 1 / np.mean(Hz)
+
+ # model function
+ # if cosine is used then tau could equal the x position of the peak
+ # if sine were to be used then tau would be the first zero crossing
+ if lock_frequency:
+ func = lambda x, A, tau: A * np.sin(2 * pi * Hz * (x - tau) + pi / 2)
+ else:
+ func = lambda x, A, Hz, tau: A * np.sin(2 * pi * Hz * (x - tau) +
+ pi / 2)
+ #func = lambda x, A, Hz, tau: A * np.cos(2 * pi * Hz * (x - tau))
+
+
+ #get peaks
+ fitted_peaks = []
+ for raw_peaks in [max_raw, min_raw]:
+ peak_data = []
+ for peak in raw_peaks:
+ index = peak[0]
+ x_data = x_axis[index - points // 2: index + points // 2 + 1]
+ y_data = y_axis[index - points // 2: index + points // 2 + 1]
+ # get a first approximation of tau (peak position in time)
+ tau = x_axis[index]
+ # get a first approximation of peak amplitude
+ A = peak[1]
+
+ # build list of approximations
+ if lock_frequency:
+ p0 = (A, tau)
+ else:
+ p0 = (A, Hz, tau)
+
+ # subtract offset from waveshape
+ y_data -= offset
+ popt, pcov = curve_fit(func, x_data, y_data, p0)
+ # retrieve tau and A i.e x and y value of peak
+ x = popt[-1]
+ y = popt[0]
+
+ # create a high resolution data set for the fitted waveform
+ x2 = np.linspace(x_data[0], x_data[-1], points * 10)
+ y2 = func(x2, *popt)
+
+ # add the offset to the results
+ y += offset
+ y2 += offset
+ y_data += offset
+
+ peak_data.append([x, y, [x2, y2]])
+
+ fitted_peaks.append(peak_data)
+
+ # structure date for output
+ max_peaks = map(lambda x: [x[0], x[1]], fitted_peaks[0])
+ max_fitted = map(lambda x: x[-1], fitted_peaks[0])
+ min_peaks = map(lambda x: [x[0], x[1]], fitted_peaks[1])
+ min_fitted = map(lambda x: x[-1], fitted_peaks[1])
+
+
+ #pylab.plot(x_axis, y_axis)
+ #pylab.hold(True)
+ #for max_p, max_f in zip(max_peaks, max_fitted):
+ # pylab.plot(max_p[0], max_p[1], 'x')
+ # pylab.plot(max_f[0], max_f[1], 'o', markersize = 2)
+ #for min_p, min_f in zip(min_peaks, min_fitted):
+ # pylab.plot(min_p[0], min_p[1], 'x')
+ # pylab.plot(min_f[0], min_f[1], 'o', markersize = 2)
+ #pylab.show()
+
+ return [max_peaks, min_peaks]
+
+
+def peakdetect_sine_locked(y_axis, x_axis, points = 9):
+ """
+ Convinience function for calling the 'peakdetect_sine' function with
+ the lock_frequency argument as True.
+
+ keyword arguments:
+ y_axis -- A list containg the signal over which to find peaks
+ x_axis -- A x-axis whose values correspond to the y_axis list and is used
+ in the return to specify the postion of the peaks.
+ points -- (optional) How many points around the peak should be used during
+ curve fitting, must be odd (default: 9)
+
+ return -- see 'peakdetect_sine'
+ """
+ return peakdetect_sine(y_axis, x_axis, points, True)
+
+
+def peakdetect_zero_crossing(y_axis, x_axis = None, window = 11):
+ """
+ Function for detecting local maximas and minmias in a signal.
+ Discovers peaks by dividing the signal into bins and retrieving the
+ maximum and minimum value of each the even and odd bins respectively.
+ Division into bins is performed by smoothing the curve and finding the
+ zero crossings.
+
+ Suitable for repeatable signals, where some noise is tolerated. Excecutes
+ faster than 'peakdetect', although this function will break if the offset
+ of the signal is too large. It should also be noted that the first and
+ last peak will probably not be found, as this function only can find peaks
+ between the first and last zero crossing.
+
+ keyword arguments:
+ y_axis -- A list containg the signal over which to find peaks
+ x_axis -- (optional) A x-axis whose values correspond to the y_axis list
+ and is used in the return to specify the postion of the peaks. If
+ omitted an index of the y_axis is used. (default: None)
+ window -- the dimension of the smoothing window; should be an odd integer
+ (default: 11)
+
+ return -- two lists [max_peaks, min_peaks] containing the positive and
+ negative peaks respectively. Each cell of the lists contains a tupple
+ of: (position, peak_value)
+ to get the average peak value do: np.mean(max_peaks, 0)[1] on the
+ results to unpack one of the lists into x, y coordinates do:
+ x, y = zip(*tab)
+ """
+ # check input data
+ x_axis, y_axis = _datacheck_peakdetect(x_axis, y_axis)
+
+ zero_indices = zero_crossings(y_axis, window = window)
+ period_lengths = np.diff(zero_indices)
+
+ bins_y = [y_axis[index:index + diff] for index, diff in
+ zip(zero_indices, period_lengths)]
+ bins_x = [x_axis[index:index + diff] for index, diff in
+ zip(zero_indices, period_lengths)]
+
+ even_bins_y = bins_y[::2]
+ odd_bins_y = bins_y[1::2]
+ even_bins_x = bins_x[::2]
+ odd_bins_x = bins_x[1::2]
+ hi_peaks_x = []
+ lo_peaks_x = []
+
+ #check if even bin contains maxima
+ if abs(even_bins_y[0].max()) > abs(even_bins_y[0].min()):
+ hi_peaks = [bin.max() for bin in even_bins_y]
+ lo_peaks = [bin.min() for bin in odd_bins_y]
+ # get x values for peak
+ for bin_x, bin_y, peak in zip(even_bins_x, even_bins_y, hi_peaks):
+ hi_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])
+ for bin_x, bin_y, peak in zip(odd_bins_x, odd_bins_y, lo_peaks):
+ lo_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])
+ else:
+ hi_peaks = [bin.max() for bin in odd_bins_y]
+ lo_peaks = [bin.min() for bin in even_bins_y]
+ # get x values for peak
+ for bin_x, bin_y, peak in zip(odd_bins_x, odd_bins_y, hi_peaks):
+ hi_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])
+ for bin_x, bin_y, peak in zip(even_bins_x, even_bins_y, lo_peaks):
+ lo_peaks_x.append(bin_x[np.where(bin_y==peak)[0][0]])
+
+ max_peaks = [[x, y] for x,y in zip(hi_peaks_x, hi_peaks)]
+ min_peaks = [[x, y] for x,y in zip(lo_peaks_x, lo_peaks)]
+
+ return [max_peaks, min_peaks]
+
+
+def _smooth(x, window_len=11, window='hanning'):
+ """
+ smooth the data using a window of the requested size.
+
+ This method is based on the convolution of a scaled window on the signal.
+ The signal is prepared by introducing reflected copies of the signal
+ (with the window size) in both ends so that transient parts are minimized
+ in the begining and end part of the output signal.
+
+ input:
+ x: the input signal
+ window_len: the dimension of the smoothing window; should be an odd
+ integer
+ window: the type of window from 'flat', 'hanning', 'hamming',
+ 'bartlett', 'blackman'
+ flat window will produce a moving average smoothing.
+
+ output:
+ the smoothed signal
+
+ example:
+
+ t = linspace(-2,2,0.1)
+ x = sin(t)+randn(len(t))*0.1
+ y = _smooth(x)
+
+ see also:
+
+ numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman,
+ numpy.convolve, scipy.signal.lfilter
+
+ TODO: the window parameter could be the window itself if a list instead of
+ a string
+ """
+ if x.ndim != 1:
+ raise ValueError, "smooth only accepts 1 dimension arrays."
+
+ if x.size < window_len:
+ raise ValueError, "Input vector needs to be bigger than window size."
+
+ if window_len<3:
+ return x
+
+ if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
+ raise(ValueError,
+ "Window is not one of '{0}', '{1}', '{2}', '{3}', '{4}'".format(
+ *('flat', 'hanning', 'hamming', 'bartlett', 'blackman')))
+
+ s = np.r_[x[window_len-1:0:-1], x, x[-1:-window_len:-1]]
+ #print(len(s))
+ if window == 'flat': #moving average
+ w = np.ones(window_len,'d')
+ else:
+ w = eval('np.' + window + '(window_len)')
+
+ y = np.convolve(w / w.sum(), s, mode = 'valid')
+ return y
+
+
+def zero_crossings(y_axis, window = 11):
+ """
+ Algorithm to find zero crossings. Smoothens the curve and finds the
+ zero-crossings by looking for a sign change.
+
+
+ keyword arguments:
+ y_axis -- A list containg the signal over which to find zero-crossings
+ window -- the dimension of the smoothing window; should be an odd integer
+ (default: 11)
+
+ return -- the index for each zero-crossing
+ """
+ # smooth the curve
+ length = len(y_axis)
+ x_axis = np.asarray(range(length), int)
+
+ # discard tail of smoothed signal
+ y_axis = _smooth(y_axis, window)[:length]
+ zero_crossings = np.where(np.diff(np.sign(y_axis)))[0]
+ indices = [x_axis[index] for index in zero_crossings]
+
+ # check if zero-crossings are valid
+ diff = np.diff(indices)
+ if diff.std() / diff.mean() > 0.2:
+ print diff.std() / diff.mean()
+ print np.diff(indices)
+ raise(ValueError,
+ "False zero-crossings found, indicates problem {0} or {1}".format(
+ "with smoothing window", "problem with offset"))
+ # check if any zero crossings were found
+ if len(zero_crossings) < 1:
+ raise(ValueError, "No zero crossings found")
+
+ return indices
+ # used this to test the fft function's sensitivity to spectral leakage
+ #return indices + np.asarray(30 * np.random.randn(len(indices)), int)
+
+############################Frequency calculation#############################
+# diff = np.diff(indices)
+# time_p_period = diff.mean()
+#
+# if diff.std() / time_p_period > 0.1:
+# raise ValueError,
+# "smoothing window too small, false zero-crossing found"
+#
+# #return frequency
+# return 1.0 / time_p_period
+##############################################################################
+
+
+
+
+
+def _test_zero():
+ _max, _min = peakdetect_zero_crossing(y,x)
+def _test():
+ _max, _min = peakdetect(y,x, delta=0.30)
+
+
+def _test_graph():
+ i = 10000
+ x = np.linspace(0,3.7*pi,i)
+ y = (0.3*np.sin(x) + np.sin(1.3 * x) + 0.9 * np.sin(4.2 * x) + 0.06 *
+ np.random.randn(i))
+ y *= -1
+ x = range(i)
+
+ _max, _min = peakdetect(y,x,750, 0.30)
+ xm = [p[0] for p in _max]
+ ym = [p[1] for p in _max]
+ xn = [p[0] for p in _min]
+ yn = [p[1] for p in _min]
+
+ plot = pylab.plot(x,y)
+ pylab.hold(True)
+ pylab.plot(xm, ym, 'r+')
+ pylab.plot(xn, yn, 'g+')
+
+ _max, _min = peak_det_bad.peakdetect(y, 0.7, x)
+ xm = [p[0] for p in _max]
+ ym = [p[1] for p in _max]
+ xn = [p[0] for p in _min]
+ yn = [p[1] for p in _min]
+ pylab.plot(xm, ym, 'y*')
+ pylab.plot(xn, yn, 'k*')
+ pylab.show()
+
+
+
+if __name__ == "__main__":
+ from math import pi
+ import pylab
+
+ i = 10000
+ x = np.linspace(0,3.7*pi,i)
+ y = (0.3*np.sin(x) + np.sin(1.3 * x) + 0.9 * np.sin(4.2 * x) + 0.06 *
+ np.random.randn(i))
+ y *= -1
+
+ #_max, _min = peakdetect(y, x, 750, 0.30)
+ _max, _min = peakdetect(y, x, lookahead=1000)
+ xm = [p[0] for p in _max]
+ ym = [p[1] for p in _max]
+ xn = [p[0] for p in _min]
+ yn = [p[1] for p in _min]
+
+ plot = pylab.plot(x, y)
+ pylab.hold(True)
+ pylab.plot(xm, ym, 'r+')
+ pylab.plot(xn, yn, 'g+')
+
+
+ pylab.show()
diff --git a/pydsf.py b/pydsf.py
new file mode 100644
index 0000000..a900857
--- /dev/null
+++ b/pydsf.py
@@ -0,0 +1,566 @@
+#! /usr/bin/env python2
+# -*- coding: utf-8 -*-
+import os
+import csv
+
+try:
+ import matplotlib as mpl
+ mpl.use('Qt4Agg')
+ import matplotlib.ticker as ticker
+ import matplotlib.pyplot as plt
+except ImportError:
+ raise ImportError('----- Matplotlib must be installed. -----')
+
+from peakdetect import *
+
+try:
+ import numpy as np
+except ImportError:
+ raise ImportError('----- NumPy must be installed. -----')
+
+try:
+ from scipy.signal import filtfilt, butter
+ from scipy import interpolate
+ from scipy import optimize
+except ImportError:
+ raise ImportError('----- SciPy must be installed. -----')
+
+class Well:
+ def __init__(self, owner):
+ self.owner = owner
+ self.name = None
+ self.raw = np.zeros(self.owner.reads, dtype=np.float)
+ self.filtered = np.zeros(self.owner.reads, dtype=np.float)
+ self.derivatives = np.zeros((4, self.owner.reads))
+ self.splines = {"raw": None,
+ "filtered": None,
+ "derivative1": None}
+ self.tm = None
+ self.tm_sd = None
+
+ def filter_raw(self):
+ """
+ Apply a filter to the raw data
+ """
+ b, a = butter(3, 0.3)
+ self.filtered = filtfilt(b, a, self.raw)
+
+ def calc_spline(self, y):
+ """
+ Calculate a spline that represents the smoothed data points
+ """
+ spline = interpolate.InterpolatedUnivariateSpline(self.owner.temprange, y)
+ return spline
+
+ def calc_derivatives(self, spline='filtered'):
+ for t in self.owner.temprange:
+ temp = self.splines[spline].derivatives(t)
+ for i in range(4):
+ self.derivatives[i, t - self.owner.t1] = temp[i]
+
+ def calc_tm(self):
+ """
+ Calculate the Tm of the well. Returns either the Tm or 'None'.
+ """
+
+ # Check if the well has already been flagged as denatured
+ if self in self.owner.denatured_wells:
+ return None # Return 'None' if true
+
+ # First assume that the well is denatured
+ self.owner.denatured_wells.append(self)
+
+ x = self.owner.temprange
+ y = self.derivatives[1]
+ lookahead = len(x)/4 # Amount of data points to look around potential peaks to verify
+ try:
+ peaks = peakdetect(y, x, lookahead) # Run peak finding algorithm
+ except:
+ return None # In case of error, return no peak
+
+ try:
+ for i in peaks[0]: # Iterate over all peaks
+ tm = i[0] # Highest peak found is the first in list
+ # Check if the peak is within cutoff range
+ if tm and tm >= self.owner.tm_cutoff_low and tm <= self.owner.tm_cutoff_high:
+ self.owner.denatured_wells.remove(self) # If everything is fine, remove the denatured flag
+ return tm # and return the Tm
+ else:
+ return None # otherwise, return no Tm
+ except:
+ return None # In case of error, return nothing
+
+ def is_denatured(self):
+ """
+ Check if the well is denatured. Returns true if the well has been already flagged as
+ denatured, no Tm was found, or if the initial signal intensity is above a user definded
+ threshold.
+ """
+ denatured = True # Assumption is that the well is denatured
+
+ if self in self.owner.denatured_wells: # check if the well is already flagged as denatured
+ return denatured # return true if it is
+
+ if self.tm and (self.tm <= self.owner.tm_cutoff_low or self.tm >= self.owner.tm_cutoff_high):
+ denatured = True
+ return denatured
+
+ for i in self.derivatives[1]: # Iterate over all points in the first derivative
+ if i > 0: # If a positive slope is found
+ denatured = False # set denatured flag to False
+
+ reads = int(round(self.owner.reads/10)) # How many values should be checked against the signal threshold:
+ # 1/10 of the total number of data point
+ read = 0 # Initialize running variable representing the current data point
+
+ if not denatured:
+ for j in self.filtered: # Iterate over the filtered data
+ if self.owner.signal_threshold: # If a signal threshold was defined
+ if j > self.owner.signal_threshold and read <= reads: # iterate over 1/10 of all data points
+ # and check for values larger than the threshold.
+ denatured = True # Set flag to True if a match is found
+ print("{}: {}".format(self.name, j))
+ return denatured # and return
+ read += 1
+
+ return denatured
+
+ def analyze(self):
+ self.filter_raw()
+ self.splines["raw"] = self.calc_spline(self.raw)
+ self.splines["filtered"] = self.calc_spline(self.filtered)
+
+ self.calc_derivatives()
+ if self.is_denatured():
+ self.owner.denatured_wells.append(self)
+
+ self.splines["derivative1"] = self.calc_spline(self.derivatives[1])
+
+ self.tm = self.calc_tm()
+ #print(self.tm)
+
+class Experiment:
+ def __init__(self, type, gui=None, files=None, replicates=None, t1=25, t2=95, dt=1, cols=12, rows=8, cutoff_low=None, cutoff_high=None, signal_threshold=None, color_range=None):
+ self.replicates = replicates
+ self.cols = cols
+ self.rows = rows
+ self.t1 = t1
+ self.t2 = t2
+ self.dt = dt
+ self.temprange = np.arange(self.t1, self.t2 + 1, self.dt, dtype=float)
+ self.reads = int(round((t2 + 1 - t1) / dt))
+ self.wellnum = self.cols * self.rows
+ self.files = files
+ self.type = type
+ self.wells = []
+ self.max_tm = None
+ self.min_tm = None
+ self.replicates = None
+ self.gui=gui
+ self.signal_threshold = signal_threshold
+ self.avg_plate = None
+ if cutoff_low:
+ self.tm_cutoff_low = cutoff_low
+ else:
+ self.tm_cutoff_low = self.t1
+ if cutoff_high:
+ self.tm_cutoff_high = cutoff_high
+ else:
+ self.tm_cutoff_high = self.t2
+ if color_range:
+ self.color_range = color_range
+ else:
+ self.color_range = None
+
+ self.plates = []
+
+ i = 1
+ for file in files:
+ plate = Plate(type=self.type, filename=file, t1=self.t1, t2=self.t2, dt=self.dt, cols=self.cols, rows=self.rows, cutoff_low=self.tm_cutoff_low, cutoff_high=self.tm_cutoff_high, signal_threshold=self.signal_threshold, color_range=self.color_range)
+ plate.id = i
+ self.plates.append(plate)
+ i += 1
+ if len(files) > 1:
+ self.avg_plate = Plate(type=self.type, filename=None, t1=self.t1, t2=self.t2, dt=self.dt, cols=self.cols, rows=self.rows, cutoff_low=self.tm_cutoff_low, cutoff_high=self.tm_cutoff_high, signal_threshold=self.signal_threshold, color_range=self.color_range)
+ self.avg_plate.id = 'average'
+
+ def analyze(self):
+ for plate in self.plates:
+ plate.analyze(gui=self.gui)
+
+ if len(self.plates) > 1:
+
+ #self.tm_replicates = np.zeros( self.wellnum, dtype=float )
+ #self.tm_replicates_sd = np.zeros( self.wellnum, dtype=float )
+
+
+ for i in range(self.wellnum):
+ tmp = []
+ for plate in self.plates:
+ tm = plate.wells[i].tm
+ self.avg_plate.wells[i].name = plate.wells[i].name
+ if plate.wells[i] not in plate.denatured_wells:
+ tmp.append(tm)
+ if len(tmp) > 0:
+ #self.avg_plate.wells[i].tm = (sum(tmp)/len(tmp))
+ self.avg_plate.wells[i].tm = np.mean(tmp)
+ self.avg_plate.wells[i].tm_sd = np.std(tmp)
+ #self.tm_replicates[i] = (sum(tmp)/len(tmp))
+ else:
+ self.avg_plate.denatured_wells.append(self.avg_plate.wells[i])
+
+
+class Plate:
+ def __init__(self, type, id=None, filename=None, replicates=None, t1=25, t2=95, dt=1, cols=12, rows=8, cutoff_low=None, cutoff_high=None, signal_threshold=None, color_range=None):
+ self.cols = cols
+ self.rows = rows
+ self.t1 = t1
+ self.t2 = t2
+ self.dt = dt
+ self.temprange = np.arange(self.t1, self.t2 + 1, self.dt, dtype=float)
+ self.reads = int(round((t2 + 1 - t1) / dt))
+ self.wellnum = self.cols * self.rows
+ self.filename = filename
+ self.type = type
+ self.wells = []
+ self.max_tm = None
+ self.min_tm = None
+ self.replicates = None
+ self.signal_threshold = signal_threshold
+ self.id = id
+ if cutoff_low:
+ self.tm_cutoff_low = cutoff_low
+ else:
+ self.tm_cutoff_low = self.t1
+ if cutoff_high:
+ self.tm_cutoff_high = cutoff_high
+ else:
+ self.tm_cutoff_high = self.t2
+ if color_range:
+ self.color_range = color_range
+ else:
+ self.color_range = None
+
+ self.denatured_wells = []
+ self.tms = []
+
+ for i in range(self.wellnum):
+ well = Well(owner = self)
+ self.wells.append(well)
+
+ #self.analyze()
+
+
+ def analytikJena(self):
+ """
+ Data processing for Analytik Jena qTower 2.0 export files
+ """
+ with open(self.filename, 'r') as f:
+ reader = csv.reader(f, delimiter=';', quoting=csv.QUOTE_NONE)
+
+ i = 0
+ for row in reader:
+ temp = np.zeros(self.reads, dtype=float)
+ for read in range(self.reads+1):
+ if read > 0:
+ try:
+ temp[read - 1] = row[read]
+ except:
+ temp[read - 1] = 0.0
+ elif read == 0:
+ self.wells[i].name = row[read]
+ self.wells[i].raw = temp
+ i += 1
+
+ def analyze(self, gui=None):
+ try:
+ # Try to access data file in the given path
+ with open(self.filename) as f: pass
+ except IOError as e:
+ # If the file is not found, or not accessible: abort
+ print('Error accessing file: {}'.format(e))
+
+
+ if self.type == 'Analytik Jena qTOWER 2.0/2.2':
+ self.analytikJena()
+ if gui:
+ update_progress_bar(gui.pb, 1)
+ else:
+ # Raise exception, if the instrument's name is unknown
+ raise NameError('Unknown instrument type: {}'.format(self.type))
+
+ for well in self.wells:
+ well.analyze()
+ if gui:
+ update_progress_bar(gui.pb, 15)
+
+ self.tms.append(well.tm)
+
+ if self.replicates:
+ if self.replicates == 'rows':
+ print "rows"
+ if self.replicates == 'cols':
+ print "cols"
+
+ self.max_tm = max(self.tms)
+ self.min_tm = min(self.tms)
+
+ def write_tm_table(self, filename):
+ with open(filename, 'w') as f:
+ f.write('#{:<4s}{:>13s}\n'.format('ID', '"Tm [°C]"'))
+ for well in self.wells:
+ if well.tm == None or well in self.denatured_wells:
+ f.write('{:<5s}{:>12s}\n'.format(well.name, 'NaN'))
+ else:
+ f.write('{:<5s}{:>12s}\n'.format(well.name, str(well.tm)))
+
+ def write_avg_tm_table(self, filename):
+ with open(filename, 'w') as f:
+ f.write('#{:<4s}{:>13s}{:>13s}\n'.format('"ID"', '"Tm [°C]"', '"SD"'))
+ for well in self.wells:
+ if well.tm == None or well in self.denatured_wells:
+ f.write('{:<5s}{:>12s}{:>12s}\n'.format(well.name, 'NaN', 'NaN'))
+ else:
+ f.write('{:<5s}{:>12s}{:>12s}\n'.format(well.name, str(well.tm), str(well.tm_sd)))
+
+ def write_raw_table(self, filename):
+ with open(filename, 'w') as f:
+ f.write('#"Raw data"\n')
+ f.write('#{:<10s}'.format('"T [°C]"'))
+ for well in self.wells:
+ f.write('{:>15s}'.format(well.name))
+ f.write('\n')
+
+ i = 0
+ for t in self.temprange:
+ f.write('{:<10s}'.format(str(t)))
+ for well in self.wells:
+ d = well.raw[i]
+ f.write('{:>-15.3f}'.format(float(np.round(d, decimals=3))))
+ f.write('\n')
+ i += 1
+
+ def write_filtered_table(self, filename):
+ with open(filename, 'w') as f:
+ f.write('#"Filtered data" \n')
+ f.write('#{:<10s}'.format('"T [°C]"'))
+ for well in self.wells:
+ f.write('{:>15s}'.format(well.name))
+ f.write('\n')
+
+ i = 0
+ for t in self.temprange:
+ f.write('{:<10s}'.format(str(t)))
+ for well in self.wells:
+ d = well.filtered[i]
+ f.write('{:>-15.3f}'.format(float(np.round(d, decimals=3))))
+ f.write('\n')
+ i += 1
+
+ def write_derivative_table(self, filename):
+ with open(filename, 'w') as f:
+ f.write('#"Derivative dI/dT"\n')
+ f.write('#{:<10s}'.format('"T [°C]"'))
+ for well in self.wells:
+ f.write('{:>15s}'.format(well.name))
+ f.write('\n')
+
+ i = 0
+ for t in self.temprange:
+ f.write('{:<10s}'.format(str(t)))
+ for well in self.wells:
+ d = well.derivatives[1][i]
+ f.write('{:>-15.3f}'.format(float(np.round(d, decimals=3))))
+ f.write('\n')
+ i += 1
+
+
+def update_progress_bar(bar, value):
+ bar.setValue(value)
+
+def plot_tm_heatmap_average(experiment, gui=None):
+ """
+ Plot Tm heatmap (Fig. 1)
+ """
+ x = 1 # Position in columns
+ y = 1 # Position in rows
+ x_values = [] # Array holding the columns
+ y_values = [] # Array holding the rows
+ c_values = [] # Array holding the color values aka Tm
+
+# c = well.tm # If not, set color to Tm
+# if c < experiment.tm_cutoff_low: # Check if Tm is lower that the cutoff
+# c = experiment.tm_cutoff_low # If it is, set color to cutoff
+# elif c > experiment.tm_cutoff_high: # Check if Tm is higher that the cutoff
+# c = experiment.tm_cutoff_high # If it is, set color to cutoff
+# else: # If the plate is denatured
+# c = experiment.tm_cutoff_low # Set its color to the low cutoff
+ for c in experiment.tm_replicates:
+
+ x_values.append(x) # Add values to the respective arrays
+ y_values.append(y)
+ c_values.append(c)
+ x += 1 # Increase column by one
+ if x > experiment.cols: # If maximum column per row is reached
+ x = 1 # reset column to one
+ y += 1 # and increase row by one
+
+ fig1 = plt.figure() # new figure
+ ax1 = fig1.add_subplot(1, 1, 1) # A single canvas
+ ax1.autoscale(tight=True) # Scale plate size
+ ax1.xaxis.set_major_locator(ticker.MaxNLocator(experiment.cols + 1)) # n columns
+ ax1.yaxis.set_major_locator(ticker.MaxNLocator(experiment.rows + 1)) # n rows
+ if experiment.color_range:
+ cax = ax1.scatter(x_values, y_values, s=300, c=c_values, marker='s', vmin=experiment.color_range[0], vmax=experiment.color_range[1]) # plot wells and color using the colormap
+ else:
+ cax = ax1.scatter(x_values, y_values, s=300, c=c_values, marker='s') # plot wells and color using the colormap
+ ax1.invert_yaxis() # invert y axis to math plate layout
+ cbar = fig1.colorbar(cax) # show colorbar
+ ax1.set_xlabel('Columns') # set axis and colorbar label
+ ax1.set_ylabel('Rows')
+ ax1.set_title('$T_m$ heatmap (average)')
+ cbar.set_label(u"Temperature [°C]")
+
+
+def plot_tm_heatmap_single(plate, gui=None):
+ """
+ Plot Tm heatmap (Fig. 1)
+ """
+ x = 1 # Position in columns
+ y = 1 # Position in rows
+ x_values = [] # Array holding the columns
+ y_values = [] # Array holding the rows
+ c_values = [] # Array holding the color values aka Tm
+ for well in plate.wells: # Iterate over all wells
+ if well not in plate.denatured_wells: # Check if well is denatured (no Tm found)
+ c = well.tm # If not, set color to Tm
+ if c < plate.tm_cutoff_low: # Check if Tm is lower that the cutoff
+ c = plate.tm_cutoff_low # If it is, set color to cutoff
+ elif c > plate.tm_cutoff_high: # Check if Tm is higher that the cutoff
+ c = plate.tm_cutoff_high # If it is, set color to cutoff
+ else: # If the plate is denatured
+ c = plate.tm_cutoff_low # Set its color to the low cutoff
+ x_values.append(x) # Add values to the respective arrays
+ y_values.append(y)
+ c_values.append(c)
+ x += 1 # Increase column by one
+ if x > plate.cols: # If maximum column per row is reached
+ x = 1 # reset column to one
+ y += 1 # and increase row by one
+
+ fig1 = plt.figure() # new figure
+ ax1 = fig1.add_subplot(1, 1, 1) # A single canvas
+ ax1.autoscale(tight=True) # Scale plate size
+ ax1.xaxis.set_major_locator(ticker.MaxNLocator(plate.cols + 1)) # n columns
+ ax1.yaxis.set_major_locator(ticker.MaxNLocator(plate.rows + 1)) # n rows
+ if plate.color_range:
+ cax = ax1.scatter(x_values, y_values, s=300, c=c_values, marker='s', vmin=plate.color_range[0], vmax=plate.color_range[1]) # plot wells and color using the colormap
+ else:
+ cax = ax1.scatter(x_values, y_values, s=300, c=c_values, marker='s') # plot wells and color using the colormap
+ ax1.invert_yaxis() # invert y axis to math plate layout
+ cbar = fig1.colorbar(cax) # show colorbar
+ ax1.set_xlabel('Columns') # set axis and colorbar label
+ ax1.set_ylabel('Rows')
+ if str(plate.id) == 'average':
+ title = '$T_m$ heatmap (average)'
+ else:
+ title = '$T_m$ heatmap (plate #{})'.format(str(plate.id))
+ ax1.set_title(title)
+ cbar.set_label(u"Temperature [°C]")
+
+# if gui:
+# update_progress_bar(gui.pb, 50)
+
+def plot_derivative(plate, gui=None):
+ """
+ Plot derivatives (Fig. 2)
+ """
+ fig2 = plt.figure() # new figure
+ fig2.suptitle('Individual Derivatives (plate #{})'.format(str(plate.id))) # set title
+
+ for plot_num in range(1, plate.wellnum + 1): # iterate over all wells
+ well = plate.wells[plot_num - 1] # get single well based on current plot number
+ ax = fig2.add_subplot(plate.rows, plate.cols, plot_num) # add new subplot
+ ax.autoscale(tight=True) # scale to data
+ ax.set_title(well.name, size='xx-small') # set title of current subplot to well identifier
+
+ if well in plate.denatured_wells:
+ ax.patch.set_facecolor('#FFD6D6')
+
+ if plot_num == plate.wellnum - plate.cols + 1: # add axis label to the subplot in the bottom left corner of the figure
+ ax.set_xlabel(u'T [°C]', size='xx-small')
+ ax.set_ylabel('dI/dT', size='xx-small')
+
+ x = plate.temprange # set values for the x axis to the given temperature range
+ y = well.derivatives[1] # grab y values from the first derivative of the well
+
+ ax.xaxis.set_major_locator(ticker.MaxNLocator(4)) # only show three tickmarks on both axes
+ ax.yaxis.set_major_locator(ticker.MaxNLocator(4))
+ if well not in plate.denatured_wells: # check if well is denatured (without determined Tm)
+ tm = well.tm # if not, grab its Tm
+ else:
+ tm = None # else set Tm to None
+ if tm:
+ ax.axvline(x=tm) # plot vertical line at the Tm
+ ax.axvspan(plate.t1, plate.tm_cutoff_low, facecolor='0.8', alpha=0.5) # shade lower cutoff area
+ ax.axvspan(plate.tm_cutoff_high, plate.t2, facecolor='0.8', alpha=0.5) # shade higher cutoff area
+ for label in ax.get_xticklabels() + ax.get_yticklabels(): # set fontsize for all tick labels to xx-small
+ label.set_fontsize('xx-small')
+
+ cax = ax.plot(x, y) # plot data to the current subplot
+# if gui:
+# update_progress_bar(gui.pb, 75)
+
+def plot_raw(plate, gui=None):
+ """
+ Plot raw data (Fig. 3)
+ """
+ fig3 = plt.figure() # new figure
+ fig3.suptitle('Raw Data (plate #{})'.format(str(plate.id))) # set title
+
+ for plot_num in range(1, plate.wellnum + 1): # iterate over all wells
+ well = plate.wells[plot_num - 1] # get single well based on current plot number
+ ax = fig3.add_subplot(plate.rows, plate.cols, plot_num) # add new subplot
+ ax.autoscale(tight=True) # scale to data
+ ax.set_title(well.name, size='xx-small') # set title of current subplot to well identifier
+
+ if well in plate.denatured_wells:
+ ax.patch.set_facecolor('#FFD6D6')
+
+ if plot_num == plate.wellnum - plate.cols + 1: # add axis label to the subplot in the bottom left corner of the figure
+ ax.set_xlabel(u'T [°C]', size='xx-small')
+ ax.set_ylabel('I', size='xx-small')
+
+ x = plate.temprange # set values for the x axis to the given temperature range
+ y = well.raw # grab y values from the raw data of the well
+
+ ax.xaxis.set_major_locator(ticker.MaxNLocator(4)) # only show three tickmarks on both axes
+ ax.yaxis.set_major_locator(ticker.MaxNLocator(4))
+ ax.axvspan(plate.t1, plate.tm_cutoff_low, facecolor='0.8', alpha=0.5) # shade lower cutoff area
+ ax.axvspan(plate.tm_cutoff_high, plate.t2, facecolor='0.8', alpha=0.5) # shade higher cutoff area
+ for label in ax.get_xticklabels() + ax.get_yticklabels(): # set fontsize for all tick labels to xx-small
+ label.set_fontsize('xx-small')
+
+ cax = ax.plot(x, y) # plot data to the current subplot
+
+# if gui:
+# update_progress_bar(gui.pb, 100)
+# gui.pb.hide()
+
+def plot(experiment, gui=None):
+ for plate in experiment.plates:
+ plot_raw(plate)
+ plot_derivative(plate)
+ plot_tm_heatmap_single(plate)
+
+ if len(experiment.plates) > 1:
+ plot_tm_heatmap_single(experiment.avg_plate)
+ #plot_tm_heatmap_average(experiment)
+
+ plt.show()
+
+
+#plate = Plate('/home/alex/Dokumente/Universitaet/Promotion/Denaturierung/26092012/26092012-MG.csv', type='analytikJena', cutoff_low=35.0, cutoff_high=60.0, signal_threshold=50000, color_range=(42, 50))
+#plot(plate)
+
+
diff --git a/ui/Ui_mainwindow.py b/ui/Ui_mainwindow.py
new file mode 100644
index 0000000..b9e23dd
--- /dev/null
+++ b/ui/Ui_mainwindow.py
@@ -0,0 +1,271 @@
+# -*- coding: utf-8 -*-
+
+# Form implementation generated from reading ui file '/home/alex/Dokumente/Programmieren/PyDSF/ui/mainwindow.ui'
+#
+# Created: Tue Oct 30 14:57:02 2012
+# by: PyQt4 UI code generator 4.9.5
+#
+# WARNING! All changes made in this file will be lost!
+
+from PyQt4 import QtCore, QtGui
+
+try:
+ _fromUtf8 = QtCore.QString.fromUtf8
+except AttributeError:
+ _fromUtf8 = lambda s: s
+
+class Ui_MainWindow(object):
+ def setupUi(self, MainWindow):
+ MainWindow.setObjectName(_fromUtf8("MainWindow"))
+ MainWindow.resize(352, 548)
+ MainWindow.setLocale(QtCore.QLocale(QtCore.QLocale.English, QtCore.QLocale.UnitedStates))
+ self.centralWidget = QtGui.QWidget(MainWindow)
+ self.centralWidget.setLocale(QtCore.QLocale(QtCore.QLocale.English, QtCore.QLocale.UnitedStates))
+ self.centralWidget.setObjectName(_fromUtf8("centralWidget"))
+ self.gridLayout_2 = QtGui.QGridLayout(self.centralWidget)
+ self.gridLayout_2.setObjectName(_fromUtf8("gridLayout_2"))
+ self.groupBox_experiment = QtGui.QGroupBox(self.centralWidget)
+ self.groupBox_experiment.setAlignment(QtCore.Qt.AlignLeading|QtCore.Qt.AlignLeft|QtCore.Qt.AlignVCenter)
+ self.groupBox_experiment.setFlat(False)
+ self.groupBox_experiment.setCheckable(False)
+ self.groupBox_experiment.setObjectName(_fromUtf8("groupBox_experiment"))
+ self.gridLayout = QtGui.QGridLayout(self.groupBox_experiment)
+ self.gridLayout.setObjectName(_fromUtf8("gridLayout"))
+ self.label_instrument = QtGui.QLabel(self.groupBox_experiment)
+ self.label_instrument.setObjectName(_fromUtf8("label_instrument"))
+ self.gridLayout.addWidget(self.label_instrument, 0, 0, 1, 1)
+ self.comboBox_instrument = QtGui.QComboBox(self.groupBox_experiment)
+ sizePolicy = QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Fixed)
+ sizePolicy.setHorizontalStretch(0)
+ sizePolicy.setVerticalStretch(0)
+ sizePolicy.setHeightForWidth(self.comboBox_instrument.sizePolicy().hasHeightForWidth())
+ self.comboBox_instrument.setSizePolicy(sizePolicy)
+ self.comboBox_instrument.setObjectName(_fromUtf8("comboBox_instrument"))
+ self.comboBox_instrument.addItem(_fromUtf8(""))
+ self.gridLayout.addWidget(self.comboBox_instrument, 0, 1, 1, 1)
+ self.groupBox_data = QtGui.QGroupBox(self.groupBox_experiment)
+ self.groupBox_data.setEnabled(True)
+ self.groupBox_data.setObjectName(_fromUtf8("groupBox_data"))
+ self.gridLayout_4 = QtGui.QGridLayout(self.groupBox_data)
+ self.gridLayout_4.setObjectName(_fromUtf8("gridLayout_4"))
+ self.listWidget_data = QtGui.QListWidget(self.groupBox_data)
+ self.listWidget_data.setAlternatingRowColors(True)
+ self.listWidget_data.setObjectName(_fromUtf8("listWidget_data"))
+ self.gridLayout_4.addWidget(self.listWidget_data, 0, 0, 1, 1)
+ self.buttonBox_open = QtGui.QDialogButtonBox(self.groupBox_data)
+ sizePolicy = QtGui.QSizePolicy(QtGui.QSizePolicy.MinimumExpanding, QtGui.QSizePolicy.Fixed)
+ sizePolicy.setHorizontalStretch(0)
+ sizePolicy.setVerticalStretch(0)
+ sizePolicy.setHeightForWidth(self.buttonBox_open.sizePolicy().hasHeightForWidth())
+ self.buttonBox_open.setSizePolicy(sizePolicy)
+ self.buttonBox_open.setLayoutDirection(QtCore.Qt.LeftToRight)
+ self.buttonBox_open.setOrientation(QtCore.Qt.Horizontal)
+ self.buttonBox_open.setStandardButtons(QtGui.QDialogButtonBox.Open)
+ self.buttonBox_open.setCenterButtons(False)
+ self.buttonBox_open.setObjectName(_fromUtf8("buttonBox_open"))
+ self.gridLayout_4.addWidget(self.buttonBox_open, 0, 1, 1, 1)
+ self.groupBox_replicates = QtGui.QGroupBox(self.groupBox_data)
+ self.groupBox_replicates.setCheckable(True)
+ self.groupBox_replicates.setChecked(False)
+ self.groupBox_replicates.setObjectName(_fromUtf8("groupBox_replicates"))
+ self.gridLayout_3 = QtGui.QGridLayout(self.groupBox_replicates)
+ self.gridLayout_3.setObjectName(_fromUtf8("gridLayout_3"))
+ self.radioButton_rep_files = QtGui.QRadioButton(self.groupBox_replicates)
+ self.radioButton_rep_files.setEnabled(False)
+ self.radioButton_rep_files.setChecked(True)
+ self.radioButton_rep_files.setObjectName(_fromUtf8("radioButton_rep_files"))
+ self.gridLayout_3.addWidget(self.radioButton_rep_files, 0, 0, 1, 1)
+ self.gridLayout_4.addWidget(self.groupBox_replicates, 1, 0, 1, 2)
+ self.gridLayout.addWidget(self.groupBox_data, 1, 0, 1, 2)
+ self.groupBox_temp = QtGui.QGroupBox(self.groupBox_experiment)
+ self.groupBox_temp.setEnabled(True)
+ self.groupBox_temp.setAutoFillBackground(False)
+ self.groupBox_temp.setCheckable(False)
+ self.groupBox_temp.setObjectName(_fromUtf8("groupBox_temp"))
+ self.formLayout = QtGui.QFormLayout(self.groupBox_temp)
+ self.formLayout.setFieldGrowthPolicy(QtGui.QFormLayout.ExpandingFieldsGrow)
+ self.formLayout.setLabelAlignment(QtCore.Qt.AlignRight|QtCore.Qt.AlignTrailing|QtCore.Qt.AlignVCenter)
+ self.formLayout.setFormAlignment(QtCore.Qt.AlignHCenter|QtCore.Qt.AlignTop)
+ self.formLayout.setObjectName(_fromUtf8("formLayout"))
+ self.label_tmin = QtGui.QLabel(self.groupBox_temp)
+ self.label_tmin.setObjectName(_fromUtf8("label_tmin"))
+ self.formLayout.setWidget(0, QtGui.QFormLayout.LabelRole, self.label_tmin)
+ self.doubleSpinBox_tmin = QtGui.QDoubleSpinBox(self.groupBox_temp)
+ self.doubleSpinBox_tmin.setAlignment(QtCore.Qt.AlignLeading|QtCore.Qt.AlignLeft|QtCore.Qt.AlignVCenter)
+ self.doubleSpinBox_tmin.setDecimals(1)
+ self.doubleSpinBox_tmin.setProperty("value", 25.0)
+ self.doubleSpinBox_tmin.setObjectName(_fromUtf8("doubleSpinBox_tmin"))
+ self.formLayout.setWidget(0, QtGui.QFormLayout.FieldRole, self.doubleSpinBox_tmin)
+ self.label_tmax = QtGui.QLabel(self.groupBox_temp)
+ self.label_tmax.setObjectName(_fromUtf8("label_tmax"))
+ self.formLayout.setWidget(1, QtGui.QFormLayout.LabelRole, self.label_tmax)
+ self.doubleSpinBox_tmax = QtGui.QDoubleSpinBox(self.groupBox_temp)
+ self.doubleSpinBox_tmax.setDecimals(1)
+ self.doubleSpinBox_tmax.setProperty("value", 95.0)
+ self.doubleSpinBox_tmax.setObjectName(_fromUtf8("doubleSpinBox_tmax"))
+ self.formLayout.setWidget(1, QtGui.QFormLayout.FieldRole, self.doubleSpinBox_tmax)
+ self.label_dt = QtGui.QLabel(self.groupBox_temp)
+ self.label_dt.setObjectName(_fromUtf8("label_dt"))
+ self.formLayout.setWidget(2, QtGui.QFormLayout.LabelRole, self.label_dt)
+ self.doubleSpinBox_dt = QtGui.QDoubleSpinBox(self.groupBox_temp)
+ self.doubleSpinBox_dt.setDecimals(1)
+ self.doubleSpinBox_dt.setProperty("value", 1.0)
+ self.doubleSpinBox_dt.setObjectName(_fromUtf8("doubleSpinBox_dt"))
+ self.formLayout.setWidget(2, QtGui.QFormLayout.FieldRole, self.doubleSpinBox_dt)
+ self.gridLayout.addWidget(self.groupBox_temp, 2, 0, 1, 1)
+ self.groupBox_cutoff = QtGui.QGroupBox(self.groupBox_experiment)
+ self.groupBox_cutoff.setEnabled(True)
+ self.groupBox_cutoff.setAlignment(QtCore.Qt.AlignLeading|QtCore.Qt.AlignLeft|QtCore.Qt.AlignVCenter)
+ self.groupBox_cutoff.setCheckable(True)
+ self.groupBox_cutoff.setChecked(False)
+ self.groupBox_cutoff.setObjectName(_fromUtf8("groupBox_cutoff"))
+ self.formLayout_2 = QtGui.QFormLayout(self.groupBox_cutoff)
+ self.formLayout_2.setFieldGrowthPolicy(QtGui.QFormLayout.ExpandingFieldsGrow)
+ self.formLayout_2.setLabelAlignment(QtCore.Qt.AlignRight|QtCore.Qt.AlignTrailing|QtCore.Qt.AlignVCenter)
+ self.formLayout_2.setFormAlignment(QtCore.Qt.AlignLeading|QtCore.Qt.AlignLeft|QtCore.Qt.AlignTop)
+ self.formLayout_2.setObjectName(_fromUtf8("formLayout_2"))
+ self.label_cutoff_high = QtGui.QLabel(self.groupBox_cutoff)
+ self.label_cutoff_high.setObjectName(_fromUtf8("label_cutoff_high"))
+ self.formLayout_2.setWidget(0, QtGui.QFormLayout.LabelRole, self.label_cutoff_high)
+ self.doubleSpinBox_upper = QtGui.QDoubleSpinBox(self.groupBox_cutoff)
+ self.doubleSpinBox_upper.setPrefix(_fromUtf8(""))
+ self.doubleSpinBox_upper.setDecimals(1)
+ self.doubleSpinBox_upper.setObjectName(_fromUtf8("doubleSpinBox_upper"))
+ self.formLayout_2.setWidget(0, QtGui.QFormLayout.FieldRole, self.doubleSpinBox_upper)
+ self.label_cutoff_low = QtGui.QLabel(self.groupBox_cutoff)
+ self.label_cutoff_low.setObjectName(_fromUtf8("label_cutoff_low"))
+ self.formLayout_2.setWidget(1, QtGui.QFormLayout.LabelRole, self.label_cutoff_low)
+ self.doubleSpinBox_lower = QtGui.QDoubleSpinBox(self.groupBox_cutoff)
+ self.doubleSpinBox_lower.setDecimals(1)
+ self.doubleSpinBox_lower.setObjectName(_fromUtf8("doubleSpinBox_lower"))
+ self.formLayout_2.setWidget(1, QtGui.QFormLayout.FieldRole, self.doubleSpinBox_lower)
+ self.gridLayout.addWidget(self.groupBox_cutoff, 2, 1, 1, 1)
+ self.groupBox_signal_threshold = QtGui.QGroupBox(self.groupBox_experiment)
+ self.groupBox_signal_threshold.setEnabled(True)
+ self.groupBox_signal_threshold.setCheckable(True)
+ self.groupBox_signal_threshold.setChecked(False)
+ self.groupBox_signal_threshold.setObjectName(_fromUtf8("groupBox_signal_threshold"))
+ self.verticalLayout = QtGui.QVBoxLayout(self.groupBox_signal_threshold)
+ self.verticalLayout.setObjectName(_fromUtf8("verticalLayout"))
+ self.spinBox_signal_threshold = QtGui.QSpinBox(self.groupBox_signal_threshold)
+ self.spinBox_signal_threshold.setMaximum(1000000)
+ self.spinBox_signal_threshold.setObjectName(_fromUtf8("spinBox_signal_threshold"))
+ self.verticalLayout.addWidget(self.spinBox_signal_threshold)
+ self.gridLayout.addWidget(self.groupBox_signal_threshold, 3, 0, 1, 1)
+ self.groupBox_cbar = QtGui.QGroupBox(self.groupBox_experiment)
+ self.groupBox_cbar.setEnabled(True)
+ self.groupBox_cbar.setCheckable(True)
+ self.groupBox_cbar.setChecked(False)
+ self.groupBox_cbar.setObjectName(_fromUtf8("groupBox_cbar"))
+ self.formLayout_4 = QtGui.QFormLayout(self.groupBox_cbar)
+ self.formLayout_4.setFieldGrowthPolicy(QtGui.QFormLayout.ExpandingFieldsGrow)
+ self.formLayout_4.setObjectName(_fromUtf8("formLayout_4"))
+ self.label_cbar_start = QtGui.QLabel(self.groupBox_cbar)
+ self.label_cbar_start.setObjectName(_fromUtf8("label_cbar_start"))
+ self.formLayout_4.setWidget(0, QtGui.QFormLayout.LabelRole, self.label_cbar_start)
+ self.doubleSpinBox_cbar_start = QtGui.QDoubleSpinBox(self.groupBox_cbar)
+ self.doubleSpinBox_cbar_start.setDecimals(1)
+ self.doubleSpinBox_cbar_start.setObjectName(_fromUtf8("doubleSpinBox_cbar_start"))
+ self.formLayout_4.setWidget(0, QtGui.QFormLayout.FieldRole, self.doubleSpinBox_cbar_start)
+ self.label_cbar_end = QtGui.QLabel(self.groupBox_cbar)
+ self.label_cbar_end.setObjectName(_fromUtf8("label_cbar_end"))
+ self.formLayout_4.setWidget(2, QtGui.QFormLayout.LabelRole, self.label_cbar_end)
+ self.doubleSpinBox_cbar_end = QtGui.QDoubleSpinBox(self.groupBox_cbar)
+ self.doubleSpinBox_cbar_end.setDecimals(1)
+ self.doubleSpinBox_cbar_end.setObjectName(_fromUtf8("doubleSpinBox_cbar_end"))
+ self.formLayout_4.setWidget(2, QtGui.QFormLayout.FieldRole, self.doubleSpinBox_cbar_end)
+ self.gridLayout.addWidget(self.groupBox_cbar, 3, 1, 1, 1)
+ self.gridLayout_2.addWidget(self.groupBox_experiment, 0, 0, 1, 1)
+ self.buttonBox_process = QtGui.QDialogButtonBox(self.centralWidget)
+ self.buttonBox_process.setStandardButtons(QtGui.QDialogButtonBox.Cancel|QtGui.QDialogButtonBox.Ok)
+ self.buttonBox_process.setObjectName(_fromUtf8("buttonBox_process"))
+ self.gridLayout_2.addWidget(self.buttonBox_process, 1, 0, 1, 1)
+ MainWindow.setCentralWidget(self.centralWidget)
+ self.menuBar = QtGui.QMenuBar(MainWindow)
+ self.menuBar.setGeometry(QtCore.QRect(0, 0, 352, 24))
+ self.menuBar.setLocale(QtCore.QLocale(QtCore.QLocale.English, QtCore.QLocale.UnitedStates))
+ self.menuBar.setObjectName(_fromUtf8("menuBar"))
+ self.menuFile = QtGui.QMenu(self.menuBar)
+ self.menuFile.setLocale(QtCore.QLocale(QtCore.QLocale.English, QtCore.QLocale.UnitedStates))
+ self.menuFile.setObjectName(_fromUtf8("menuFile"))
+ self.menuHelp = QtGui.QMenu(self.menuBar)
+ self.menuHelp.setLocale(QtCore.QLocale(QtCore.QLocale.English, QtCore.QLocale.UnitedStates))
+ self.menuHelp.setObjectName(_fromUtf8("menuHelp"))
+ MainWindow.setMenuBar(self.menuBar)
+ self.statusBar = QtGui.QStatusBar(MainWindow)
+ self.statusBar.setObjectName(_fromUtf8("statusBar"))
+ MainWindow.setStatusBar(self.statusBar)
+ self.actionQuit = QtGui.QAction(MainWindow)
+ self.actionQuit.setObjectName(_fromUtf8("actionQuit"))
+ self.actionAbout = QtGui.QAction(MainWindow)
+ self.actionAbout.setObjectName(_fromUtf8("actionAbout"))
+ self.actionAbout_Qt = QtGui.QAction(MainWindow)
+ icon = QtGui.QIcon()
+ icon.addPixmap(QtGui.QPixmap(_fromUtf8(":/qtlogo.svg")), QtGui.QIcon.Normal, QtGui.QIcon.Off)
+ self.actionAbout_Qt.setIcon(icon)
+ self.actionAbout_Qt.setObjectName(_fromUtf8("actionAbout_Qt"))
+ self.menuFile.addAction(self.actionQuit)
+ self.menuHelp.addAction(self.actionAbout)
+ self.menuHelp.addAction(self.actionAbout_Qt)
+ self.menuBar.addAction(self.menuFile.menuAction())
+ self.menuBar.addAction(self.menuHelp.menuAction())
+ self.label_instrument.setBuddy(self.comboBox_instrument)
+ self.label_tmin.setBuddy(self.doubleSpinBox_tmin)
+ self.label_tmax.setBuddy(self.doubleSpinBox_tmax)
+ self.label_dt.setBuddy(self.doubleSpinBox_dt)
+ self.label_cutoff_high.setBuddy(self.doubleSpinBox_upper)
+ self.label_cutoff_low.setBuddy(self.doubleSpinBox_lower)
+ self.label_cbar_start.setBuddy(self.doubleSpinBox_cbar_start)
+ self.label_cbar_end.setBuddy(self.doubleSpinBox_cbar_end)
+
+ self.retranslateUi(MainWindow)
+ QtCore.QMetaObject.connectSlotsByName(MainWindow)
+
+ def retranslateUi(self, MainWindow):
+ MainWindow.setWindowTitle(QtGui.QApplication.translate("MainWindow", "PyDSF", None, QtGui.QApplication.UnicodeUTF8))
+ self.groupBox_experiment.setTitle(QtGui.QApplication.translate("MainWindow", "Experimental Setup", None, QtGui.QApplication.UnicodeUTF8))
+ self.label_instrument.setText(QtGui.QApplication.translate("MainWindow", "Instrument", None, QtGui.QApplication.UnicodeUTF8))
+ self.comboBox_instrument.setItemText(0, QtGui.QApplication.translate("MainWindow", "Analytik Jena qTOWER 2.0/2.2", None, QtGui.QApplication.UnicodeUTF8))
+ self.groupBox_data.setToolTip(QtGui.QApplication.translate("MainWindow", "Add data files to the experiment. If multiple files are loaded, they are treated as replicates.
", None, QtGui.QApplication.UnicodeUTF8))
+ self.groupBox_data.setTitle(QtGui.QApplication.translate("MainWindow", "Data File", None, QtGui.QApplication.UnicodeUTF8))
+ self.groupBox_replicates.setTitle(QtGui.QApplication.translate("MainWindow", "Replicates", None, QtGui.QApplication.UnicodeUTF8))
+ self.radioButton_rep_files.setText(QtGui.QApplication.translate("MainWindow", "Files", None, QtGui.QApplication.UnicodeUTF8))
+ self.groupBox_temp.setToolTip(QtGui.QApplication.translate("MainWindow", "