qdsfit/QDS.py
Markus Rosenstihl 805a726346 Fit limits are shown and labeles
Size of PeakBox elements are now OK
Conductivity moved to data
2013-07-10 18:36:07 +02:00

649 lines
24 KiB
Python
Executable File

#!/usr/bin/env python
# -*- encoding: utf-8 -*-
import os
import sys
import re
import signal
from PyQt4.QtCore import *
from PyQt4.QtGui import *
#from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas
from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar
#import matplotlib.pyplot as plt
import matplotlib.colors
import numpy as N
import scipy.odr as O
import scipy.optimize as opt
import MatplotlibWidget
import QDSMain
import PeakWidget
from data import Data, Conductivity, conductivity
#import yaff
def sigint_handler(*args):
"""Handler for the SIGINT signal (CTRL + C).
"""
sys.stderr.write('\r')
if QMessageBox.question(None, '', "Are you sure you want to quit?",
QMessageBox.Yes | QMessageBox.No,
QMessageBox.No) == QMessageBox.Yes:
QApplication.quit()
def id_to_color(id):
"""
"""
colors = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
conv = matplotlib.colors.ColorConverter()
return conv.to_rgb(colors[id % len(colors)])
def tau_peak(f, a, b):
tau = (N.sin(N.pi * a / 2. / (b + 1)) / N.sin(N.pi * a * b / 2. / (b + 1))) ** (1 / a)
tau /= 2 * N.pi * f
return tau
def hn(p, nu):
delta_eps, tau, a, b = p
om = 2 * N.pi * nu
Phi = N.arctan((om * tau) ** a * N.sin(N.pi * a / 2.) / (1. + (om * tau) ** a * N.cos(N.pi * a / 2.)))
e_loss = delta_eps * (1 + 2 * (om * tau) ** a * N.cos(N.pi * a / 2.) + (om * tau) ** (2. * a) ) ** (
-b / 2.) * N.sin(b * Phi)
e_stor = delta_eps * (1 + 2 * (om * tau) ** a * N.cos(N.pi * a / 2.) + (om * tau) ** (2. * a) ) ** (
-b / 2.) * N.cos(b * Phi)
return 2 * e_loss
def multi_hn(p, nu):
conductivity = p[1]
cond_beta = p[2]
om = 2 * N.pi * nu
e_loss = conductivity / om ** cond_beta
e_loss += p[0]
#for key, igroup in groupby(p[3:], lambda x: x//4):
for i in xrange(len(p[3:]) / 4):
delta_eps, tau, a, b = p[3 + i * 4:3 + (i + 1) * 4]
#delta_eps, tau, a, b = list(igroup)
#print delta_eps,tau,a,b
#a = 0.5 *(1 + N.tanh(a))
#b = 0.5 *(1 + N.tanh(b))
Phi = N.arctan((om * tau) ** a * N.sin(N.pi * a / 2.) / (1. + (om * tau) ** a * N.cos(N.pi * a / 2.)))
e_loss += 2 * delta_eps * (1 + 2 * (om * tau) ** a * N.cos(N.pi * a / 2.) + (om * tau) ** (2. * a) ) ** (
-b / 2.) * N.sin(b * Phi)
#e_stor = delta_eps * (1+ 2*(om*tau)**a * N.cos(N.pi*a/2.) + (om*tau)**(2.*a) )**(-b/2.)*N.cos(b*Phi)
return e_loss
def mini_func(p, x, y):
res = y - multi_hn(p, x)
# apply weights
res /= 1 / y
return N.sqrt(N.dot(res, res))
def fit_odr(x, y, p0, fixed):
dat = O.Data(x, y, 1.0 / y)
mod = O.Model(multi_hn)
odr = O.ODR(dat, mod, p0, ifixx=(0,), ifixb=fixed, maxit=2000)
odr.run()
return odr.output.beta
def fit_lbfgsb(x, y, p0, fixed):
# TODO fixed parameters…
bounds = [(0, None), (0, 1)]
for i in xrange(len(p0[3:]) / 4):
bounds.append((1e-4, 1e12)) # delta_eps
bounds.append((1e-12, 1e3)) # tau
bounds.append((0.1, 1)) # a
bounds.append((0.1, 1)) # b
x, f_minvalue, info_dict = opt.fmin_l_bfgs_b(mini_func, p0,
fprime=None,
args=(x, y),
approx_grad=True,
bounds=bounds,
iprint=0,
maxfun=4000)
if info_dict['warnflag'] != 0:
print info_dict["task"]
return x
def fit_anneal(x, y, p0, fixed):
bounds = [(0, 1e14), (0, 1)]
for i in xrange(len(p0[2:]) / 4):
bounds.append((1e-4, 1e12)) # delta_eps
bounds.append((1e-12, 1e3)) # tau
bounds.append((0.1, 1)) # a
bounds.append((0.1, 1)) # b
ret = opt.anneal(mini_func, p0,
args=(x, y),
maxeval=20000,
maxiter=30000,
lower=[b[0] for b in bounds],
upper=[b[1] for b in bounds],
dwell=100,
full_output=1)
#pmin, func_min, final_Temp, cooling_iters,accepted_tests, retval
#retval : int
#Flag indicating stopping condition::
# 0 : Points no longer changing
# 1 : Cooled to final temperature
# 2 : Maximum function evaluations
# 3 : Maximum cooling iterations reached
# 4 : Maximum accepted query locations reached
# 5 : Final point not the minimum amongst encountered points
print "Stop reason", ret
return ret[0]
class Peak(QObject):
changedData = pyqtSignal()
def __init__(self, id=None, mpl=None):
QObject.__init__(self)
super(Peak, self).__init__()
self.color = id_to_color(id)
self.widget = PeakWidget.PeakWidget()
self.widget.setId(id)
self.widget.setColor(map(int, [255 * i for i in self.color]))
self.widget.changedTable.connect(self.updatePeak)
self.mpl = mpl
self.mpl_line = None
def getParameter(self):
p = self.widget.peakParameter()
return p
def getFixed(self):
p = self.widget.fixedParameter()
return p
def setParameter(self, delta_eps=None, tau=None, a=None, b=None):
self.widget.updateTable(delta_eps, tau, a, b)
self.updatePeak()
def updatePeak(self):
# get current axis limits
x_min, x_max = self.mpl.canvas.axes.get_xlim()
y_min, y_max = self.mpl.canvas.axes.get_ylim()
nu = N.logspace(N.log10(x_min), N.log10(x_max), 2048)
y = hn(self.getParameter(), nu)
# clip data to axes limits
mask = (y < y_max) & (y > y_min)
y = y[mask]
nu = nu[mask]
if self.mpl_line == None:
self.mpl_line, = self.mpl.canvas.axes.loglog(nu, y, '--', label="Peak %i" % (self.widget.id),
animated=True) # peak
self.mpl_line.set_color(self.color)
else:
self.mpl_line.set_xdata(nu)
self.mpl_line.set_ydata(y)
self.changedData.emit()
class AppWindow(QMainWindow):
def __init__(self, parent=None):
super(AppWindow, self).__init__(parent)
self.ui = QDSMain.Ui_MainWindow()
self.ui.setupUi(self)
self.picked_artist = None
self.data = None
self.Conductivity = None
self._lines = dict()
self.peakId = 0
self.peakBoxes = {}
## menubar
fileMenu = self.menuBar().addMenu("File")
openFile = QAction("&Open", self)
openFile.setShortcut(QKeySequence.Open)
openFile.triggered.connect(self.openFile)
fileMenu.addAction(openFile)
# fitting methods
fitMenu = self.menuBar().addMenu("Standard Fits")
# lm
fit_lmAction = QAction("&Levenberg-Marquardt", self)
fit_lmAction.setShortcut(QKeySequence("Ctrl+F"))
fitMenu.addAction(fit_lmAction)
# lbfgsb
fit_lbfgsbAction = QAction("&L-BFGS-B", self)
fitMenu.addAction(fit_lbfgsbAction)
# Simulated Annealing
fit_annealAction = QAction("&Simulated Annealing", self)
fitMenu.addAction(fit_annealAction)
# YAFF
yaffMenu = self.menuBar().addMenu("YAFF")
start_yaff = QAction("&Startparam.", self)
yaffMenu.addAction(start_yaff)
fit_yaff = QAction("&Fit", self)
yaffMenu.addAction(fit_yaff)
self.signalMapper = QSignalMapper(self)
for i, fit_action in enumerate([fit_lmAction, fit_lbfgsbAction, fit_annealAction
]):
self.signalMapper.setMapping(fit_action, i)
fit_action.triggered.connect(self.signalMapper.map)
self.signalMapper.mapped.connect(self.fitData)
# save fitted values
self.ui.actionSave_FitResult.triggered.connect(self.saveFit)
# the plot area, a matplotlib widget
self.mplWidget = PlotWidget(self.ui.mplWidget)
self.mplWidget.canvas.draw()
self.mplWidget.updateGeometry()
# what to do with CIDs?
self.cid = []
self.cid.append(self.mplWidget.canvas.mpl_connect("button_press_event", self.mpl_button_press))
self.cid.append(self.mplWidget.canvas.mpl_connect("pick_event", self.mpl_button_pick))
self.cid.append(self.mplWidget.canvas.mpl_connect("button_release_event", self.mpl_button_release))
#self.cid.append(self.mplWidget.canvas.mpl_connect("resize_event", self.myresizeEvent))
self.mplWidget.toolbar.spanSelectedTrigger.connect(self.set_fit_xlimits)
def resizeEvent(self, evt):
"""resizing the main window executes this method
TODO: FIX RESIZEING keeping mainwindow maximised for now
:param evt:
"""
self.mplWidget.canvas.draw()
self.mplWidget._bg_cache = self.mplWidget.canvas.copy_from_bbox(self.mplWidget.canvas.axes.bbox)
for line in self.mplWidget.canvas.axes.get_lines():
line.set_animated(False)
self.mplWidget.canvas.draw()
for line in self.mplWidget.canvas.axes.get_lines():
line.set_animated(True)
def saveFit(self):
"""
Saving fit parameters to fitresults.log
including temperature
"""
if not os.path.exists("fitresults.log"):
f = open("fitresults.log", "w")
# write header
f.write('# T ')
if self.Conductivity != None:
f.write("%8s %8s %8s " % ("e_s", "sig", "pow_sig"))
for i, pb in enumerate(self.peakBoxes):
enum_peak = ("e_inf_%i" % i, "tau_%i" % i, "alpha_%i" % i, "beta_%i" % i)
f.write("%8s %8s %8s %8s " % enum_peak)
f.write("high_lim lower_lim")
f.write('\n')
f.flush()
else:
f = open("fitresults.log", "a")
#f.write("%3.2f "%(self.data.meta["T"]))
pars = list(self.fitresult)
pars.insert(0, self.data.meta["T"])
N.savetxt(f, N.array([pars, ]), fmt="%8.5g", delimiter=" ")
f.close()
def set_fit_xlimits(self, xmin, xmax):
self.data.fit_limits = (xmin, xmax, None, None)
self.updatePlot()
def mpl_button_release(self, event):
ax = self.mplWidget.canvas.axes
if self.picked_artist:
if not event.inaxes: # moved outside the plot, add back to original position
ax.add_artist(self.picked_artist)
else: # we move one of the three points determinig the peak
self.picked_artist.set_xdata(event.xdata)
self.picked_artist.set_ydata(event.ydata)
ax.add_artist(self.picked_artist)
for peak in self.peakBoxes.keys():
peak.updatePeak()
self.picked_artist = None
self.mplWidget.canvas.draw_idle()
def mpl_button_pick(self, event):
self.picked_artist = event.artist
event.artist.remove()
self.mplWidget.canvas.draw_idle()
def addCond(self, pos):
if self.Conductivity != None:
return
self.statusBar().showMessage("Click on graph")
self.Conductivity = Conductivity(mpl=self.mplWidget)
self.Conductivity.changedData.connect(self.updatePlot)
self.Conductivity.setParameter(0, 1 / (pos[0] / pos[1] / 2 / N.pi), 1.0)
self.ui.scrollAreaWidgetContents.layout().addWidget(self.Conductivity.widget)
self.Conductivity.widget.ui.removeButton.clicked.connect(self.delCond)
def delCond(self):
self.cond_param = None
self.cond = None
self.Conductivity.mpl_line.remove()
self.Conductivity.mpl_line_static.remove()
del self.Conductivity
self.Conductivity = None
self.updatePlot()
def addPeak(self, pos):
self.peakId += 1
self.statusBar().showMessage("Select Peak Position")
peak = Peak(id=self.peakId, mpl=self.mplWidget)
# connect to delPeak
peak.widget.ui.removeButton.clicked.connect(self.delPeak)
peak.setParameter(delta_eps=pos[1], tau=1 / (2. * N.pi * pos[0]), a=1, b=1)
peak.changedData.connect(self.updatePlot)
self.peakBoxes[peak] = None
for pb in self.peakBoxes.keys():
self.ui.scrollAreaWidgetContents.layout().addWidget(pb.widget)
self.updatePlot()
def delPeak(self):
deletePeaks = []
for i in xrange(self.ui.scrollAreaWidgetContents.layout().count()):
print i
for i, peak in enumerate(self.peakBoxes.keys()):
if peak.widget.isHidden():
peak.mpl_line.remove()
deletePeaks.append(peak)
for peak in deletePeaks:
self.peakBoxes.pop(peak)
self.updatePlot()
def fitData(self, method):
if self.Conductivity != None:
start_parameter = list(self.Conductivity.getParameter())
fixed_params = [i for i in self.Conductivity.getFixed()]
else:
start_parameter = [0, 0, 1]
fixed_params = [0, 0, 0]
for pb in self.peakBoxes.keys():
[start_parameter.append(i) for i in pb.getParameter()]
[fixed_params.append(i) for i in pb.getFixed()]
fit_methods = [fit_odr, fit_lbfgsb, fit_anneal]
print "StartParameter", start_parameter
print "FixedParameter", fixed_params
print "Limits (xmin, xmax, ymin, ymax)", self.data.fit_limits
_freq, _fit = self.data.get_data()
result = fit_methods[method](_freq, _fit.imag, start_parameter, fixed_params)
self.fitresult = result
for i, pb in enumerate(self.peakBoxes.keys()):
delta_eps, tau, a, b = result[3 + i * 4:3 + (i + 1) * 4]
pb.setParameter(delta_eps, tau, a, b)
e_static, sigma, sigma_N = result[:3]
if self.Conductivity != None:
self.Conductivity.setParameter(e_static, sigma, sigma_N)
#print "*** FIT RESULTS ***"
#print u"\u03c3"
#print u"\u0394\u03b5"
self.updatePlot()
def mpl_button_press(self, event):
"""
Handles the clicks on the matplotlib figure canvas
"""
if self.ui.actionAdd_Peak.isChecked() and event.inaxes:
x, y = event.xdata, event.ydata
self.addPeak((x, y))
self.ui.actionAdd_Peak.setChecked(False)
self.statusBar().clear()
if self.ui.actionAdd_Cond.isChecked() and event.inaxes:
x, y = event.xdata, event.ydata
self.addCond((x, y))
self.ui.actionAdd_Cond.setChecked(False)
self.statusBar().clear()
def openFile(self):
ax = self.mplWidget.canvas.axes
path = unicode(QFileDialog.getOpenFileName(self, "Open file"))
# TODO anaylize file (LF,MF, HF) and act accordingly
data = N.loadtxt(path, skiprows=4)
numpat = re.compile('\d+\.\d+')
try:
Temp = None
for line in open(path).readlines():
if re.search("Fixed", line) or re.search("Temp", line):
print "Found line with Fixed or Temp"
Temp = float(re.search(numpat, line).group())
print "Temperature found in file:", Temp
break
if Temp == None: raise ValueError
except:
Temp = QInputDialog.getDouble(self, "No temperature found in data set", "Temperature/K:", value=0.00)[0]
# mask the data to values > 0 (loglog plot)
mask = (data[:, 1] > 0) & (data[:, 2] > 0) #& (data[:,2]>1e-3) & (data[:,0] > 1e-2)
_freq = data[mask, 0]
_die_stor = data[mask, 1]
_die_loss = data[mask, 2]
# clear the figure
ax.clear()
#if self.data != None:
# self.data.remove_curves()
self.data = Data(_freq, _die_stor, _die_loss)
self.data.meta["T"] = Temp
self.data.data_curve, = ax.loglog(self.data.frequency,
self.data.epsilon.imag,
'b.',
markersize=4,
label="Data",
animated=True)
ax.set_xlabel("Frequency/Hz", fontsize=16)
#ax.set_ylabel(u'\u03B5"', fontsize=16)
ax.set_ylabel(u'e"', fontsize=16)
ax.autoscale(True)
ax.set_xlim(_freq.min()/3, _freq.max()*3)
self.legend = ax.legend(title="T=%.2f" % Temp)
for line in ax.get_lines():
line.set_animated(False)
ax.grid()
self.mplWidget.canvas.draw()
# weird behaviour need to draw all first, then set_animated the redraw
for line in ax.get_lines():
line.set_animated(True)
self.legend.set_animated(True)
ax.autoscale(False)
self.mplWidget._bg_cache = self.mplWidget.canvas.copy_from_bbox(ax.bbox)
self.updatePlot()
def updatePlot(self):
ax = self.mplWidget.canvas.axes
nu = self.data.frequency
fit = N.zeros(len(nu))
for peak in self.peakBoxes.keys():
params = peak.getParameter()
fit += hn(params, nu)
if self.Conductivity != None:
print "Cond. given"
params = self.Conductivity.getParameter()[1:]
fit += conductivity(params, nu)
fit += self.Conductivity.getParameter()[0] # eps static
# clip data to axes limits
y_min, y_max = ax.get_ylim()
mask = (fit < y_max) & (fit > y_min)
#mask = N.ones(len(fit), dtype="bool")
if self.data.fitted_curve == None:
self.data.fitted_curve, = ax.loglog(nu[mask], fit[mask],
'k-',
alpha=0.5,
label="Sum",
animated=True)
else:
self.data.fitted_curve.set_xdata(nu[mask])
self.data.fitted_curve.set_ydata(fit[mask])
if self.mplWidget._axvlims != []:
[axv.remove() for axv in self.mplWidget._axvlims]
self.mplWidget._axvlims = []
if self.mplWidget._axvname != []:
[axvname.remove() for axvname in self.mplWidget._axvname]
self.mplWidget._axvname = []
for xlim in self.data.fit_limits[:2]:
self.mplWidget._axvlims.append(ax.axvline(xlim, color="k", ls="--", alpha=0.5, animated=True))
self.mplWidget._axvname.append(ax.text(xlim, y_min*3. , "%.3g"%xlim,
horizontalalignment='center',
verticalalignment='center',
animated=True) )
self.mplWidget.canvas.restore_region(self.mplWidget._bg_cache)
self.legend = ax.legend(title="T=%.2f" % (self.data.meta["T"]))
self.legend.set_animated(True)
for animated_artist in ax.findobj(match=lambda x: x.get_animated()):
#print "updatePlot: animated artist:",animated_artist
ax.draw_artist(animated_artist)
self.mplWidget.canvas.blit(ax.bbox)
class PlotWidget(QWidget):
def __init__(self, parent=None):
QWidget.__init__(self)
super(PlotWidget, self).__init__(parent)
self.mplwidget = MatplotlibWidget.MatplotlibWidget(hold=True,
xlim=(1e-2, 1e7),
xscale='log',
yscale='log')
self.canvas = self.mplwidget.figure.canvas # shortcut
self.canvas.axes.grid(True)
#self.bbox_size = self.canvas.axes.bbox.size
self.toolbar = CustomToolbar(self.canvas, self.mplwidget, parent)
self.toolbar.setToolButtonStyle(Qt.ToolButtonTextUnderIcon)
layout = QVBoxLayout(parent)
#self.mplwidget.setLayout(layout)
layout.addWidget(self.canvas)
layout.addWidget(self.mplwidget)
layout.addWidget(self.toolbar)
self._bg_cache = None
self._axvlims = []
self._axvname = []
class CustomToolbar(NavigationToolbar):
# our spanChanged signal
spanSelectedTrigger = pyqtSignal(float, float, name='spanChanged')
def __init__(self, plotCanvas, plotWidget, parent=None):
NavigationToolbar.__init__(self, plotCanvas, plotWidget, coordinates=True)
self.canvas = plotCanvas
# Span select Button
#self.span_button = QAction(QIcon("border-1d-right-icon.png" ), "Span", self)
self.span_button = QAction(QIcon(QPixmap(":/icons/fit_limits.png")), "Fit Limits", self)
self.span_button.setCheckable(True)
self.cids = []
self.cids.append(self.canvas.mpl_connect('button_press_event', self.press))
self.cids.append(self.canvas.mpl_connect('motion_notify_event', self.onmove))
self.cids.append(self.canvas.mpl_connect('button_release_event', self.release))
self.cids.append(self.canvas.mpl_connect('draw_event', self.update_background))
# act.setCheckable(True)
# add actions before the coordinates widget
self.insertAction(self.actions()[-1], self.span_button)
self.buttons = {}
self._pressed = False
self.background = None
self.span = None
self.istart = 0
self.iend = 0
self.xstart = 0
self.xend = 0
def set_span(self, x1, x2):
#trans = blended_transform_factory(self.axes.transData, self.axes.transAxes)
cur = self.span.get_xy()
cur[0, 0] = x1
cur[1, 0] = x1
cur[2, 0] = x2
cur[3, 0] = x2
self.span.set_xy(cur)
def ignore(self, event):
# 'return ``True`` if *event* should be ignored'
return event.inaxes != self.canvas.axes or event.button != 1
def update_background(self, event):
#if self.canvas.axes is None:
# raise SyntaxError,"Need an axes reference!"
self.background = self.canvas.copy_from_bbox(self.canvas.axes.bbox)
def press(self, event):
if self.span_button.isChecked():
if self.background is None:
self.update_background()
else:
self.canvas.restore_region(self.background)
self.xstart = event.xdata
self.istart = event.x
if self.span is None:
self.span = self.canvas.axes.axvspan(event.xdata, event.xdata, alpha=0.10, color="k", animated=False)
else:
self.set_span(event.xdata, event.xdata)
self._pressed = True
def onmove(self, event):
if self.span_button.isChecked() and self._pressed and not self.ignore(event):
self.set_span(self.xstart, event.xdata)
self.update()
def update(self):
self.canvas.restore_region(self.background)
self.canvas.axes.draw_artist(self.span)
for line in self.canvas.axes.get_lines():
self.canvas.axes.draw_artist(line)
self.canvas.blit(self.canvas.axes.bbox)
def release(self, event):
self.span_button.setChecked(False)
self.xend = event.xdata
self.iend = event.x
if self.iend < self.istart:
self.iend, self.istart = self.istart, self.iend
print "released", self.xstart, self.xend
if self._pressed:
if self.ignore(event):
self.istart = 0
self.spanSelectedTrigger.emit(self.xstart, self.xend)
self._pressed = False
if __name__ == '__main__':
signal.signal(signal.SIGINT, sigint_handler)
app = QApplication(sys.argv)
timer = QTimer()
timer.start(1000) # Check every second for Strg-c on Cmd line
timer.timeout.connect(lambda: None)
main = AppWindow()
main.showMaximized()
main.raise_()
sys.exit(app.exec_())