qdsfit/QDS.py

550 lines
20 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 *
import matplotlib
from mathlib import fit_anneal, fit_lbfgsb, fit_odr_cmplx, FunctionRegister, FitRoutine
matplotlib.use('agg')
from matplotlib import pyplot
from matplotlib.colors import hex2color
#matplotlib.rc_file("default.mplrc")
import numpy as np
import QDSMain
from data import Data, Conductivity, Peak, PowerComplex, Static, YAFF
import pyqtgraph as pg
#import yaff
from CustomWidgets import ParameterWidget, YaffWidget
USE_CROSSH=False
class AppWindow(QMainWindow):
def __init__(self, parent=None):
super(AppWindow, self).__init__(parent)
self.ui = QDSMain.Ui_MainWindow()
self.ui.setupUi(self)
self._last_written_header = None
actions = {
self.ui.actionAdd_Peak:self.addPeak,
self.ui.actionAdd_Cond:self.addCond,
self.ui.actionAdd_PowerLaw:self.addPowerComplex,
self.ui.actionAdd_Eps_Infty:self.addEpsInfty,
self.ui.actionYAFF:self.addYaff,
}
self.myActionGroup = QActionGroup(self)
for a in actions.keys(): self.myActionGroup.addAction(a)
self.function_registry = FunctionRegister()
self.peakId = 0
self.parameterWidget = ParameterWidget()
self.ui.dockWidget_3.setWidget(self.parameterWidget)
## menubar
fileMenu = self.menuBar().addMenu("File")
openFile = QAction("&Open", self)
openFile.setShortcut(QKeySequence.Open)
openFile.triggered.connect(self.openFile)
fileMenu.addAction(openFile)
saveFile = QAction("&Save Fit Result", self)
saveFile.setShortcut(QKeySequence.Save)
saveFile.triggered.connect(self.saveFitResult)
fileMenu.addAction(saveFile)
# 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._fit_thread = QThread()
# self._fit_method = FitRoutine()
# self._fit_method.moveToThread(self._fit_thread)
# self._fit_method.finished.connect(self._fit_thread.quit)
# self._fit_queue = QEventLoop()
# self._fit_timer = QTimer()
# self._fit_timer.setSingleShot(True)
# self._fit_timer.timeout.connect(self._fit_queue.quit)
# self._fit_method.finished.connect(self._fit_queue.quit)
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_start)
# save fitted values
#self.ui.actionSave_FitResult.triggered.connect(self.saveFitResult)# replaced by menu
self.data = Data()
self.fit_boundary = pg.LinearRegionItem(brush=QColor(0,127,254,15))
self.ui.pgPlotWidget_imag.addItem(self.data.data_curve_imag)
self.ui.pgPlotWidget_real.addItem(self.data.data_curve_real)
self.ui.pgPlotWidget_imag.addItem(self.data.fitted_curve_imag)
self.ui.pgPlotWidget_real.addItem(self.data.fitted_curve_real)
# cross hair
if USE_CROSSH:
#self.poslabel = pg.LabelItem(justify='right')
self.horizontalLine = pg.InfiniteLine(angle=0)
self.verticalLine = pg.InfiniteLine(angle=90)
#self.ui.pgPlotWidget_imag.plotItem.addItem(self.poslabel)
self.ui.pgPlotWidget_imag.addItem(self.horizontalLine)
self.ui.pgPlotWidget_imag.addItem(self.verticalLine)
self.ui.pgPlotWidget_imag.addItem(self.fit_boundary)
for pltwidgt in (self.ui.pgPlotWidget_real, self.ui.pgPlotWidget_imag):
pltwidgt.setLogMode(x=True, y=True)
pltwidgt.showGrid(x=True, y=True)
pltwidgt.disableAutoRange()
pltwidgt.setLabel("bottom", "Frequency", units="Hz")
self.ui.pgPlotWidget_imag.setLabel("left", u'Dielectric loss ε"' , units="Debye")
self.ui.pgPlotWidget_real.setLabel("left", u"Dielectric loss ε' ", units="Debye")
sc_real = self.ui.pgPlotWidget_real.scene()
sc_real.sigMouseClicked.connect(self.mousePress)
sc_real.sigMouseMoved.connect(self.updateCrosshair)
sc_imag = self.ui.pgPlotWidget_imag.scene()
sc_imag.sigMouseClicked.connect(self.mousePress)
sc_imag.sigMouseMoved.connect(self.updateCrosshair)
def updateCrosshair(self,evt):
vb_real = self.ui.pgPlotWidget_real.getPlotItem().vb
vb_imag = self.ui.pgPlotWidget_imag.getPlotItem().vb
if self.ui.pgPlotWidget_imag.underMouse():
pos = vb_imag.mapSceneToView(evt)
elif self.ui.pgPlotWidget_real.underMouse():
pos = vb_real.mapSceneToView(evt)
else:
pos = QPointF(0.0, 0.0)
self.last_pos = pos
if USE_CROSSH:
self.horizontalLine.setBounds([self.data.frequency.min(), self.data.frequency.max()])
self.horizontalLine.setPos(pos.y())
self.verticalLine.setPos(pos.x())
#self.poslabel.setText("<span style='font-size: 12pt'>f=%0.1f Hz <span style='color: red'>e=%0.1f</span>"\
# % (pos.x(),pos.y()))
def mousePress(self, evt):
data_pos = self.last_pos
mouse_in_imag = self.ui.pgPlotWidget_imag.underMouse()
mouse_in_real = self.ui.pgPlotWidget_real.underMouse()
msgBox = QMessageBox()
if self.ui.actionAdd_Peak.isChecked():
if mouse_in_imag:
self.addPeak(data_pos)
self.ui.actionAdd_Peak.setChecked(False)
else:
msgBox.setText("Click in imaginary part")
msgBox.exec_()
if self.ui.actionAdd_Cond.isChecked():
if mouse_in_imag:
self.addCond(data_pos)
self.ui.actionAdd_Cond.setChecked(False)
else:
msgBox.setText("Click in imaginary part")
msgBox.exec_()
if self.ui.actionYAFF.isChecked():
if mouse_in_imag:
self.addYaff(data_pos)
self.ui.actionYAFF.setChecked(False)
else:
msgBox.setText("Click in imaginary part")
msgBox.exec_()
if self.ui.actionAdd_PowerLaw.isChecked():
if mouse_in_imag:
self.addPowerComplex(data_pos)
self.ui.actionAdd_PowerLaw.setChecked(False)
else:
msgBox.setText("Click in imaginary part")
msgBox.exec_()
if self.ui.actionAdd_Eps_Infty.isChecked():
if mouse_in_real:
self.addEpsInfty(data_pos)
self.ui.actionAdd_Eps_Infty.setChecked(False)
else:
msgBox.setText("Click in real part")
msgBox.exec_()
def saveFitResult(self):
"""
Saving fit parameters to fitresults.log
including temperature
"""
self.saveFitFigure()
if not os.path.exists("fitresults.log"):
f = open("fitresults.log", "w")
else:
f = open("fitresults.log", "a")
# prepare header
header="# T "
pars = []
for fcn in self.function_registry.get_registered_functions():
for name in fcn.widget.names: # variable names
header += "{n:11}{n_sd:11}".format(n=name, n_sd=name+"_sd")
for i,par in enumerate(fcn._beta): # params # TODO: ughh
if fcn._selector_mask is not None:
if fcn._selector_mask[i]:
pars.extend([par])
pars.extend([fcn._sd_beta[i]])
else:
pars.extend([par])
pars.extend([fcn._sd_beta[i]])
header += "%-11s%-11s\n"%("fit_xlow","fit_xhigh")
# write new header if fit model changed
if self._last_written_header != header:
f.write(header)
f.flush()
self._last_written_header = header
else:
pass
pars.insert(0, self.data.meta["T"])
pars.append(self.data.fit_limits[0])
pars.append(self.data.fit_limits[1])
pars = np.array([pars])
np.savetxt(f, pars, fmt = '%-10.3e', delimiter=" ")
f.close()
def saveFitFigure(self):
fig = pyplot.figure(figsize=(3.54*1.4, 2.75*1.4))
font = {'family' : 'sans serif',
'weight' : 'normal',
'size' : 9}
matplotlib.rc('font', **font)
pyplot.grid(linestyle="solid",alpha=0.3, color="0.5")
pyplot.loglog(self.data.frequency, self.data.epsilon.imag, 'bo', markersize=4, label="Data")
pyplot.loglog(self.data.frequency_fit, self.data.epsilon_fit.imag, 'r-', lw=1.2, label="Fit")
for fcn in self.function_registry.get_registered_functions():
f,eps = fcn.get_data()
label = fcn.id_string
color = hex2color(str(fcn.color.name()))
pyplot.loglog(f,eps[1], ls=":", color=color, lw=1, dashes=(1,1), label=label)
for i in (0,1): pyplot.axvline(x=self.data.fit_limits[i], color='b', ls="-", lw=0.5)
pyplot.legend(title = "T=%.1f K"%(self.data.meta["T"]))
pyplot.xlabel('f/Hz')
pyplot.ylabel(u'ε"')
pyplot.ylim(self.data.epsilon.imag.min(), self.data.epsilon.imag.max() )
#pyplot.savefig(os.path.splitext(self.filepath)[0]+".png")
pyplot.savefig(os.path.splitext(self.filepath)[0]+".pdf")
fig.clear()
def addYaff(self, pos):
_yaff = YAFF(plt_real=self.ui.pgPlotWidget_real,
plt_imag=self.ui.pgPlotWidget_imag,
limits=self.data.fit_limits)
_yaff.blockSignals(True)
_yaff.changedData.connect(self.updatePlot)
_yaff.removeObj.connect(self.delParamterObject)
gg_y = 10**pos.y()*2
gg_x = 1/(10**pos.x()*2*np.pi)
yaff_par = [ gg_y, gg_x , 20.0, 1.0, 0.5, gg_x/100, 1.0, 1.0]
_yaff.setParameter(beta=yaff_par)
self.parameterWidget.add(_yaff.widget)
self.function_registry.register_function(_yaff)
self.updatePlot()
_yaff.blockSignals(False)
def addCond(self, pos):
_conductivity = Conductivity(plt_real=self.ui.pgPlotWidget_real,
plt_imag=self.ui.pgPlotWidget_imag,
limits=self.data.fit_limits)
_conductivity.blockSignals(True)
_conductivity.changedData.connect(self.updatePlot)
_conductivity.removeObj.connect(self.delParamterObject)
cond_par = [0.0, 10**(pos.y() + pos.x())*2*np.pi , 1.0]
_conductivity.setParameter(beta=cond_par)
self.parameterWidget.add(_conductivity.widget)
self.function_registry.register_function(_conductivity) ##todo
self.updatePlot()
_conductivity.blockSignals(False)
def addPowerComplex(self, pos):
_power_complex = PowerComplex(plt_imag=self.ui.pgPlotWidget_imag,
plt_real=self.ui.pgPlotWidget_real,
limits=self.data.fit_limits)
_power_complex.changedData.connect(self.updatePlot)
_power_complex.removeObj.connect(self.delParamterObject)
cond_par = [10**(pos.y() + pos.x())*2*np.pi , 1.0]
_power_complex.setParameter(beta=cond_par)
self.parameterWidget.add(_power_complex.widget)
self.function_registry.register_function(_power_complex)
self.updatePlot()
def addEpsInfty(self, pos):
_eps_infty = Static(plt_imag=self.ui.pgPlotWidget_imag,
plt_real=self.ui.pgPlotWidget_real,
limits=self.data.fit_limits)
_eps_infty.changedData.connect(self.updatePlot)
_eps_infty.removeObj.connect(self.delParamterObject)
cond_par = [10**pos.y()]
_eps_infty.setParameter(beta=cond_par)
self.parameterWidget.add(_eps_infty.widget)
self.function_registry.register_function(_eps_infty)
self.updatePlot()
def delParamterObject(self, obj):
print "unregister",obj
self.function_registry.unregister_function(obj)
self.updatePlot()
def addPeak(self, pos):
id_list = [ key.id_num for key in
self.function_registry.get_registered_functions().keys()
if key.id_string == 'hn']
self.peakId = 1
while self.peakId in id_list:
self.peakId += 1
_peak = Peak(id_num=self.peakId,
plt_real=self.ui.pgPlotWidget_real,
plt_imag=self.ui.pgPlotWidget_imag,
limits=self.data.fit_limits)
_peak.changedData.connect(self.updatePlot)
_peak.removeObj.connect(self.delParamterObject)
new_peak = [2*10**pos.y(), 1 / (2*np.pi*10**pos.x()), 1, 1]
_peak.setParameter(beta = new_peak)
self.function_registry.register_function(_peak)
self.parameterWidget.add(_peak.widget)
self.updatePlot()
def delPeak(self):
for i, peak in enumerate(self.peakBoxes.keys()):
if peak.widget.isHidden():
self.ui.pgPlotWidget_imag.removeItem(peak.mpl_line)
self.parameterWidget.vlayout.removeWidget(peak.widget)
self.function_registry.unregister_function(peak)
self.peakBoxes.pop(peak)
self.parameterWidget.vlayout.update()
self.updatePlot()
@pyqtSlot(np.ndarray)
def test(self, t=None):
print "test",t
def fitData_start(self, method):
fit_methods = [fit_odr_cmplx, fit_lbfgsb, fit_anneal]
# build function list
p0,funcs,fixed_params = [],[],[]
for fcn in self.function_registry.get_registered_functions():
p0.extend(fcn.getParameter())
funcs.append(fcn)
fixed_params.extend(fcn.getFixed())
fcn.clearData()
_freq, _fit = self.data.get_data()
# Prepare
if not self._fit_thread.isRunning():
self._fit_method = FitRoutine()
self._fit_method.moveToThread(self._fit_thread)
self._fit_method.finished_fit.connect(self.fitData_update)
self._fit_method.fit_odr_cmplx(_freq, _fit, p0, fixed_params, funcs)
self._fit_method.data_ready.connect(self.updateIntermediatePlot)
#for fun in funcs:
# self._fit_method.data_ready.connect(fun.updateData)
self._fit_thread.started.connect(self._fit_method.fit)
self._fit_thread.finished.connect(self._fit_method.deleteLater)
self._fit_thread.start()
self.ui.statusbar.showMessage("Fitting ...")
else:
print "fit is still running"
pass
def fitData_update(self):
print "before thread.quit"
print "running thread?:",self._fit_thread.isRunning()
self._fit_thread.quit()
print "after thread.quit"
print "running thread?:",self._fit_thread.isRunning()
odr_result = self._fit_method.result()
print odr_result
# build function list
p0,funcs,fixed_params = [],[],[]
for fcn in self.function_registry.get_registered_functions():
p0.extend(fcn.getParameter())
funcs.append(fcn)
fixed_params.extend(fcn.getFixed())
#print "Set fit data"
self.data.set_fit(odr_result.beta, funcs)
self.ui.statusbar.showMessage(" ".join(odr_result.stopreason))
i = 0
for fcn in self.function_registry.get_registered_functions():
num_p = len(fcn.getParameter())
beta = odr_result.beta[i:num_p+i]
sd_beta = odr_result.sd_beta[i:num_p+i]
fcn.setParameter(beta, sd_beta)
i += num_p
#self.updatePlot()
def openFile(self):
#path = unicode(QFileDialog.getOpenFileName(self, "Open file"))
path = "MCM42PG0_199.96K.dat"
self.filepath=path
# TODO analyize file (LF,MF, HF) and act accordingly
data = np.loadtxt(path, skiprows=4)
self.setWindowTitle(os.path.basename(path))
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]
self.data.set_data(_freq, _die_stor, _die_loss)
self.data.meta["T"] = Temp
self.fit_boundary.setRegion([np.log10(_freq.min()), np.log10(_freq.max())])
self.ui.pgPlotWidget_imag.enableAutoRange()
self.ui.pgPlotWidget_real.enableAutoRange()
self.updatePlot()
self.ui.pgPlotWidget_imag.disableAutoRange()
self.ui.pgPlotWidget_real.disableAutoRange()
def updatePlot(self):
log10fmin, log10fmax = self.fit_boundary.getRegion()
self.data.set_fit_xlimits(10**log10fmin, 10**log10fmax)
p0,funcs = [],[]
for fcn in self.function_registry.get_registered_functions():
p0.extend(fcn.getParameter())
funcs.append(fcn)
# calculate parametrized curve
self.data.set_fit(p0, funcs)
# replot data and fit, TODO: replot only if data changed
self.data.data_curve_real.setData(self.data.frequency, self.data.epsilon.real)
self.data.data_curve_imag.setData(self.data.frequency, self.data.epsilon.imag)
if len(funcs)> 0:
self.data.fitted_curve_real.setData(self.data.frequency_fit, self.data.epsilon_fit.real)
self.data.fitted_curve_imag.setData(self.data.frequency_fit, self.data.epsilon_fit.imag)
else:
self.data.fitted_curve_real.setData([np.nan],[np.nan])
self.data.fitted_curve_imag.setData([np.nan],[np.nan])
def updateIntermediatePlot(self, freq, intermediate_data):
self.data.fitted_curve_real.setData(freq, intermediate_data[0])
self.data.fitted_curve_imag.setData(freq, intermediate_data[1])
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.Yes) == QMessageBox.Yes:
QApplication.quit()
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_())