From 8fc316f2ff5612e51fb55e4b0fdd72a2e3f0f7d6 Mon Sep 17 00:00:00 2001 From: Markus Rosenstihl Date: Wed, 17 Jul 2013 13:00:24 +0200 Subject: [PATCH] * moved math functions (fit, hn, etc.) to mathlib.py * fitresults are stored in a better format --- MatplotlibWidget.py | 28 ++++- QDS.py | 258 +++++--------------------------------------- QDSToolbar.py | 98 +++++++++++++++++ mathlib.py | 113 +++++++++++++++++++ 4 files changed, 265 insertions(+), 232 deletions(-) create mode 100644 QDSToolbar.py create mode 100644 mathlib.py diff --git a/MatplotlibWidget.py b/MatplotlibWidget.py index 4f1d285..3d866fe 100644 --- a/MatplotlibWidget.py +++ b/MatplotlibWidget.py @@ -15,11 +15,12 @@ This software is licensed under the terms of the MIT License Derived from 'embedding_in_pyqt4.py': Copyright © 2005 Florent Rougon, 2006 Darren Dale """ +from QDSToolbar import CustomToolbar __version__ = "1.0.0" -from PyQt4.QtGui import QSizePolicy -from PyQt4.QtCore import QSize +from PyQt4.QtGui import QSizePolicy, QWidget, QVBoxLayout +from PyQt4.QtCore import QSize, Qt from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as Canvas from matplotlib.figure import Figure @@ -123,3 +124,26 @@ if __name__ == '__main__': win = ApplicationWindow() win.show() sys.exit(app.exec_()) + + +class PlotWidget(QWidget): + def __init__(self, parent=None): + QWidget.__init__(self) + super(PlotWidget, self).__init__(parent) + self.mplwidget = 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 = [] \ No newline at end of file diff --git a/QDS.py b/QDS.py index f9e1b07..d8931f2 100755 --- a/QDS.py +++ b/QDS.py @@ -8,19 +8,13 @@ import signal from PyQt4.QtCore import * from PyQt4.QtGui import * -import matplotlib -import matplotlib.colors -#from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas -from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar -#import matplotlib.pyplot as plt +from mathlib import fit_anneal, fit_lbfgsb, fit_odr, hn, id_to_color +from matplotlibWidget import PlotWidget -matplotlib.rc_file("default.mplrc") +#matplotlib.rc_file("default.mplrc") 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 @@ -38,13 +32,7 @@ def sigint_handler(*args): 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): @@ -53,104 +41,25 @@ def tau_peak(f, a, b): 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): @@ -250,7 +159,7 @@ class AppWindow(QMainWindow): self.signalMapper.mapped.connect(self.fitData) # save fitted values - self.ui.actionSave_FitResult.triggered.connect(self.saveFit) + self.ui.actionSave_FitResult.triggered.connect(self.saveFitResult) # the plot area, a matplotlib widget self.mplWidget = PlotWidget(self.ui.mplWidget) self.mplWidget.canvas.draw() @@ -278,29 +187,32 @@ class AppWindow(QMainWindow): line.set_animated(True) - def saveFit(self): + def saveFitResult(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"])) + # write header + f.write('# T ') + parfmt = "%.2f" # T formatting + # if self.Conductivity != None: pass# always true + f.write("%8s %8s %8s " % ("e_s", "sig", "pow_sig")) + parfmt += " %.3g %.3g %.2f " # conductivity formatting + 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) + parfmt += " %.3g %.3g %.2f %.2f" # peak formatting + f.write("high_lim lower_lim") # TODO: store limits + f.write('\n') + f.flush() + #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=" ") + N.savetxt(f, N.array([pars, ]), fmt=parfmt, delimiter=" ") f.close() def set_fit_xlimits(self, xmin, xmax): @@ -532,120 +444,6 @@ class AppWindow(QMainWindow): 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) diff --git a/QDSToolbar.py b/QDSToolbar.py new file mode 100644 index 0000000..b7ce845 --- /dev/null +++ b/QDSToolbar.py @@ -0,0 +1,98 @@ +from PyQt4.QtCore import pyqtSignal +from PyQt4.QtGui import QAction, QIcon, QPixmap +from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar, NavigationToolbar2QTAgg + +__author__ = 'markusro' + + +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): + """Overrides method of NavigationToolbar. + Allows fast drawing of the span selector with blitting + """ + 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 \ No newline at end of file diff --git a/mathlib.py b/mathlib.py new file mode 100644 index 0000000..42fcdc7 --- /dev/null +++ b/mathlib.py @@ -0,0 +1,113 @@ +# -*- encoding: utf-8 -*- +import matplotlib +import numpy as N +from scipy import optimize as opt, optimize, odr +#from QDS import mini_func, multi_hn + +__author__ = 'markusro' + + +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] + + +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_odr(x, y, p0, fixed): + dat = odr.Data(x, y, 1.0 / y) + mod = odr.Model(multi_hn) + fit = odr.ODR(dat, mod, p0, ifixx=(0,), ifixb=fixed, maxit=2000) + fit.run() + return fit.output.beta + + +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 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 mini_func(p, x, y): + res = y - multi_hn(p, x) + # apply weights + res /= 1 / y + return N.sqrt(N.dot(res, res)) + + +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 \ No newline at end of file