#!/usr/bin/env python # -*- encoding: utf-8 -*- import os,sys import re import signal from PyQt4.QtCore import * from PyQt4.QtGui import * from PyQt4.uic 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 import ConductivityWidget from data import Data #import yaff def sigint_handler(*args): """Handler for the SIGINT signal.""" 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 conductivity(p, nu): c = p[0]/(2*N.pi*nu)**p[1] return c def hnfit(delta_eps, tau, points): x = [point[0] for point in points] y = [point[1] for point in points] dat = O.Data(x,y) mod = O.Model(hn) odr = O.ODR(dat,mod, [delta_eps, tau, 1, 1], ifixb=(0,0,1,1), ifixx=(0,), maxit=1000) odr.run() odr.output.pprint() return odr.output.beta 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 Conductivity(QObject): changedData = pyqtSignal() def __init__(self, mpl=None): QObject.__init__(self) super(Conductivity, self) self.widget = ConductivityWidget.ConductivityWidget() self.widget.changedTable.connect(self.updateData) self.mpl_line = None self.mpl_line_static = None self.mpl = mpl def getParameter(self): p = self.widget.getTable() return p def getFixed(self): p = self.widget.fixedParameter() return p def setParameter(self, eps_static=None, sigma=None, sigma_N=None): self.widget.updateTable(eps_static, sigma, sigma_N) self.updateData() def updateData(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), 1024) eps_static, sigma, sigma_N = self.getParameter() y = conductivity([sigma, sigma_N],nu) y_static = N.ones(len(nu))*eps_static # clip data to axes limits mask_static = (y_static < y_max) & (y_static > y_min) # clip data to axes limits mask = (y < y_max) & (y > y_min) #mask = mask_static = N.ones(1024, dtype='bool') if self.mpl_line == None: self.mpl_line, = self.mpl.canvas.axes.loglog(nu[mask],y[mask],'k--',label="Cond.", animated=True) # peak else: self.mpl_line.set_xdata(nu[mask]) self.mpl_line.set_ydata(y[mask]) if self.mpl_line_static == None: self.mpl_line_static, = self.mpl.canvas.axes.loglog(nu[mask_static], y_static[mask_static], 'k:', label=r"$\epsilon_S$", animated=True) # peak else: self.mpl_line_static.set_xdata(nu[mask_static]) self.mpl_line_static.set_ydata(y_static[mask_static]) self.changedData.emit() 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("Fit") # lm fit_lmAction = QAction("&Levenberg-Marquardt",self) 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) 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) # he 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.onclick) ) self.cid.append( self.mplWidget.canvas.mpl_connect("pick_event", self.pick) ) self.cid.append( self.mplWidget.canvas.mpl_connect("button_release_event", self.mpl_button_release) ) self.mplWidget.toolbar.spanSelectedTrigger.connect(self.set_fit_xlimits) def resizeEvent(self,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): if not os.path.exists("fitresults.log"): f = open("fitresults.log","w") # write header f.write('# T ') if self.Conductivity != None: f.write("e_s sig pow_sig ") for i,pb in enumerate(self.peakBoxes): f.write("e_inf_%i tau_%i alpha_%i beta_%i"%(i,i,i,i)) 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"] ) print pars N.savetxt(f,N.array([pars,]),fmt="%.5g", delimiter=" ") f.close() def set_fit_xlimits(self,xmin,xmax): self.data.fit_limits = (xmin, xmax, None,None) def mpl_button_release(self,event): if self.picked_artist: if not event.inaxes: # moved outside the plot, add back to original position self.mplWidget.canvas.axes.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) self.mplWidget.canvas.axes.add_artist(self.picked_artist) for peak in self.peakBoxes.keys(): peak.updatePeak() self.picked_artist = None self.mplWidget.canvas.draw_idle() def 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 onclick(self, event): """ Handles the clicks on the matplotlib figure canvas """ if self.ui.actionAdd_Peak.isChecked(): 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(): x,y = event.xdata,event.ydata self.addCond((x,y)) self.ui.actionAdd_Cond.setChecked(False) self.statusBar().clear() def openFile(self): 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: for line in open(path).readlines(): if re.search("Fixed", line) or re.search("Temp", line): Temp = float(re.search(numpat, line).group()) print "Temperature found in file:",Temp break else: print "No Temperature found in file" #Temp = QInputDialog.getDouble(self, "No temperature found in data set","Temperature/K:", value=Temp)[0] except: Temp = QInputDialog.getDouble(self, "No temperature found in data set","Temperature/K:", value=0.0)[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 self.mplWidget.canvas.axes.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, = self.mplWidget.canvas.axes.loglog(self.data.frequency, self.data.epsilon.imag, 'b.', markersize=4, label="Data", animated=True) self.mplWidget.canvas.axes.set_xlim(_freq.min(), _freq.max()) self.mplWidget.canvas.axes.set_xlabel("frequency/Hz", fontsize=16) self.mplWidget.canvas.axes.set_ylabel(u'\u03B5"', fontsize=16) self.mplWidget.canvas.axes.autoscale(True) self.legend = self.mplWidget.canvas.axes.legend(title="T=%.2f"%Temp) for line in self.mplWidget.canvas.axes.get_lines(): line.set_animated(False) self.mplWidget.canvas.axes.grid() self.mplWidget.canvas.draw() for line in self.mplWidget.canvas.axes.get_lines(): line.set_animated(True) self.legend.set_animated(True) self.mplWidget.canvas.axes.autoscale(False) self.mplWidget._bg_cache = self.mplWidget.canvas.copy_from_bbox(self.mplWidget.canvas.axes.bbox) self.updatePlot() def updatePlot(self): 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 = self.mplWidget.canvas.axes.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, = self.mplWidget.canvas.axes.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]) # update lines self.mplWidget.canvas.restore_region(self.mplWidget._bg_cache) self.legend = self.mplWidget.canvas.axes.legend(title = "T=%.2f"%(self.data.meta["T"]) ) self.legend.set_animated(True) for animated_artist in self.mplWidget.canvas.axes.findobj(match=lambda x: x.get_animated()): #print "updatePlot: animated artist:",animated_artist self.mplWidget.canvas.axes.draw_artist(animated_artist) self.mplWidget.canvas.blit( self.mplWidget.canvas.axes.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 print "PlotWidget",self.bbox_size #self.toolbar = NavigationToolbar(self.mplwidget, parent) self.toolbar = CustomToolbar(self.canvas, self.mplwidget, parent) layout = QVBoxLayout(parent) layout.addWidget(self.canvas) layout.addWidget(self.mplwidget) layout.addWidget(self.toolbar) self._bg_cache = None 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")), "Span", 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_())