diff --git a/QDS.py b/QDS.py index 69ab7f4..9e81a07 100755 --- a/QDS.py +++ b/QDS.py @@ -9,7 +9,7 @@ from PyQt4.QtCore import * from PyQt4.QtGui import * import matplotlib -from mathlib import fit_anneal, fit_lbfgsb, fit_odr, hn +from mathlib import fit_anneal, fit_lbfgsb, fit_odr, hn, FitFunctionCreator, fit_odr_cmplx matplotlib.use('agg') @@ -90,7 +90,8 @@ class AppWindow(QMainWindow): self.data = Data() self.fit_boundary = pg.LinearRegionItem(brush=QColor(254,254,254,10)) - self.ui.graphicsView.addItem(self.data.data_curve) + self.ui.graphicsView.addItem(self.data.data_curve_imag) + self.ui.graphicsView.addItem(self.data.data_curve_real) self.ui.graphicsView.addItem(self.data.fitted_curve) self.ui.graphicsView.addItem(self.fit_boundary) self.ui.graphicsView.setLogMode(x=True, y=True) @@ -124,58 +125,59 @@ class AppWindow(QMainWindow): else: f = open("fitresults.log", "a") # write header - f.write('# T ') - parfmt = "%.2f" # T formatting + f.write("#%7s"%('T')) + parfmt = "%8.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 + f.write("%9s%9s%9s " % ("e_s", "sig", "pow_sig")) + parfmt += "%9.3g%9.3g%9.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("%9s%9s%9s%9s " % enum_peak) + print enum_peak + parfmt += "%9.3g%9.3g%9.2f%9.2f" # peak formatting + f.write("fit_xlow fit_xhigh") # TODO: store limits + parfmt += "%9.3g%9.3g" f.write('\n') f.flush() #f.write("%3.2f "%(self.data.meta["T"])) pars = list(self.fitresult) pars.insert(0, self.data.meta["T"]) + pars.append(self.data.fit_limits[0]) + pars.append(self.data.fit_limits[1]) N.savetxt(f, N.array([pars, ]), fmt=parfmt, delimiter=" ") f.close() def saveFitFigure(self): - fig = pyplot.Figure(figsize=(3.54, 2.75)) + fig = pyplot.figure(figsize=(3.54, 2.75)) font = {'family' : 'sans serif', 'weight' : 'normal', - 'size' : 16} + 'size' : 8} matplotlib.rc('font', **font) - print self.data.epsilon_fit.shape, type(self.data.epsilon_fit) - pyplot.loglog(self.data.frequency, self.data.epsilon.imag, 'bo', markersize=4, label="Data") - pyplot.loglog(self.data.frequency, self.data.epsilon_fit, 'r-', lw=2, label="Fit") + + pyplot.loglog(self.data.frequency, self.data.epsilon.imag, 'bo', markersize=3, label="Data") + pyplot.loglog(self.data.frequency, self.data.epsilon_fit, 'r-', lw=1, label="Fit") for i,peak in enumerate(self.peakBoxes): f,eps = peak.get_data() color = hex2color(str(peak.get_color().name())) - pyplot.loglog(f,eps, ls="--", color=color , lw=1.5, label="Peak %i"%i) + pyplot.loglog(f,eps, ls="--", color=color , lw=0.75, label="Peak %i"%i) if self.Conductivity != None: f,eps = self.Conductivity.get_conductivity() color = hex2color(str(self.Conductivity.get_color().name())) - pyplot.loglog(f,eps, ls="-.", color=color, lw=1.5, label="Cond.") + pyplot.loglog(f,eps, ls="-.", color=color, lw=0.75, label="Cond.") f,eps = self.Conductivity.get_epsilon_static() - pyplot.loglog(f,eps, ls=":", color=color, lw=1.5, label=r'$\epsilon_0$') + pyplot.loglog(f,eps, ls=":", color=color, lw=0.75, label=r'$\epsilon_0$') - pyplot.legend(title="T=%.1f K"%(self.data.meta["T"]) ) + for i in (0,1): pyplot.axvline(x=self.data.fit_limits[i], color='g', ls="--") + pyplot.legend(title = "T=%.1f K"%(self.data.meta["T"])) pyplot.grid() pyplot.xlabel('f/Hz') pyplot.ylabel('eps"') - pyplot.savefig("test.png") + pyplot.savefig(os.path.splitext(self.filepath)[0]+".png") fig.clear() - def set_fit_xlimits(self, xmin, xmax): - self.data.fit_limits = (xmin, xmax, None, None) - self.updatePlot() - def addCond(self, pos): if self.Conductivity != None: @@ -239,18 +241,29 @@ class AppWindow(QMainWindow): [start_parameter.append(i) for i in pb.getParameter()] [fixed_params.append(i) for i in pb.getFixed()] + log10fmin, log10fmax = self.fit_boundary.getRegion() - xmin,xmax,ymin,ymax = self.data.fit_limits - self.data.fit_limits = [10**log10fmin, 10**log10fmax,ymin,ymax] + self.data.set_fit_xlimits(10**log10fmin, 10**log10fmax) 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) + + # check new method + if 1: + funcs = ["static","conductivity"] if self.Conductivity != None else [] + for pb in self.peakBoxes.keys(): + funcs.append("hn") + newres = fit_odr_cmplx(_freq, _fit, start_parameter, fixed_params, funcs) + + print newres + + self.fitresult = result for i, pb in enumerate(self.peakBoxes.keys()): - delta_eps, tau, a, b = result[3 + i * 4:3 + (i + 1) * 4] + 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: @@ -262,6 +275,7 @@ class AppWindow(QMainWindow): def openFile(self): path = unicode(QFileDialog.getOpenFileName(self, "Open file")) + self.filepath=path #path = "MCM42PG0_199.96K.dat" # TODO anaylize file (LF,MF, HF) and act accordingly data = N.loadtxt(path, skiprows=4) @@ -304,7 +318,8 @@ class AppWindow(QMainWindow): fit += self.Conductivity.getParameter()[0] # eps static self.data.epsilon_fit = fit[:] - self.data.data_curve.setData(self.data.frequency, self.data.epsilon.imag) + self.data.data_curve_imag.setData(self.data.frequency, self.data.epsilon.imag) + self.data.data_curve_imag.setData(self.data.frequency, self.data.epsilon.real) if len(self.peakBoxes) > 0 and self.Conductivity != None: self.data.fitted_curve.setData(nu, fit) diff --git a/data.py b/data.py index ae515c2..14e4c23 100644 --- a/data.py +++ b/data.py @@ -15,12 +15,14 @@ class Data: self.epsilon_fit = die_real*0 + 1j * die_imag*0 myPen = pg.mkPen(width=3, color=(255,255,127)) - self.data_curve = pg.PlotDataItem(x=[N.nan], y=[N.nan],pen=QColor(0,0,0,0), symbol='o', + self.data_curve_imag = pg.PlotDataItem(x=[N.nan], y=[N.nan],pen=QColor(0,0,0,0), symbol='o', + symbolBrush=(255,127,0,127)) + self.data_curve_real = pg.PlotDataItem(x=[N.nan], y=[N.nan],pen=QColor(0,0,0,0), symbol='s', symbolBrush=(255,127,0,127)) self.fitted_curve = pg.PlotDataItem(N.array([N.nan]), N.array([N.nan]), pen=myPen) self.length = len(frequency) self.meta = dict() - self.fit_limits = (frequency.min(), frequency.max(), die_imag.min(), die_imag.max()) + self.fit_limits = [frequency.min(), frequency.max(), die_imag.min(), die_imag.max()] def __del__(self): #self.remove_curves() @@ -31,8 +33,18 @@ class Data: self.frequency = f self.epsilon = e_real + 1j*e_imag self.epsilon_fit = 0*e_real + 1j*e_imag*0 - self.fit_limits = (f.min(), f.max(), e_imag.min(), e_imag.max()) - self.data_curve.setData(f,e_imag) + self.fit_limits = [f.min(), f.max(), e_imag.min(), e_imag.max()] + self.data_curve_imag.setData(f,e_imag) + self.data_curve_real.setData(f,e_real) + + def set_fit_xlimits(self, xmin, xmax): + self.fit_limits[0] = xmin + self.fit_limits[1] = xmax + + def set_fit_ylimits(self, ymin, ymax): + self.fit_limits[2] = ymin + self.fit_limits[3] = ymax + def get_data(self): """ @@ -40,7 +52,8 @@ class Data: """ mask = N.ones(len(self.frequency), dtype='bool') mask = (self.frequency > self.fit_limits[0]) & (self.frequency < self.fit_limits[1]) - mask &= (self.epsilon.imag > self.fit_limits[2]) & (self.epsilon.imag < self.fit_limits[1]) + mask &= (self.epsilon.imag > self.fit_limits[2]) & (self.epsilon.imag < self.fit_limits[3]) + mask &= (self.epsilon.real > self.fit_limits[2]) & (self.epsilon.real < self.fit_limits[3]) return self.frequency[mask], self.epsilon[mask] def remove_curves(self): diff --git a/mathlib.py b/mathlib.py index c6818b4..f768939 100644 --- a/mathlib.py +++ b/mathlib.py @@ -124,4 +124,80 @@ def multi_hn(p, nu): 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 \ No newline at end of file + return tau + + +### define funcs here +class Functions: + def __init__(self): + self.list = { + "hn":(self.hn_cmplx,4), + "conductivity":(self.cond_cmplx,1), + "power":(self.power_cmplx,2), + "static":(self.static_cmplx,1), + } + + def hn_cmplx(self, p, x): + om = 2*N.pi*x + hn = om*1j + eps,t,a,b = p + hn = eps/(1+(1j*om*t)**a)**b + cplx = N.array([hn.real, -hn.imag]) + return cplx + + def cond_cmplx(self, p, x): + om = 2*N.pi*x + sgma = p[0] + cond = sgma/(1j*om) + cplx = N.array([cond.real, -cond.imag]) + return cplx + + def power_cmplx(self, p, x): + om = 2*N.pi*x + sgma,n = p + power = sgma/(om*1j)**n + cplx = N.array([power.real, -power.imag]) + return cplx + + def static_cmplx(self, p, x): + eps_inf = p[0] + static = N.ones((2, len(x)))*eps_inf + static[:,1] *= 0 # set imag part zero + #cplx = N.array([static.real, static.imag]) + return static + + def get(self,name): + return self.list[name] + + +class FitFunctionCreator: + def __init__(self): + self.data = None + self.functions = Functions() + + def fitfcn(self, p0, x, *funcs): + self.data = N.zeros( x.shape ) + ndx = 0 + for fn in funcs: # loop over functions and add the results up + f,num_p = self.functions.get(fn) + p = p0[ndx:ndx+num_p] + self.data += f(p, x[0]) # fit functions take only 1-dim x + ndx += num_p + return self.data + +def fit_odr_cmplx(x, y, p0, fixed, fcns): + f = FitFunctionCreator() + if x.ndim < 2: + x = N.resize(x, (2,x.size)) + if N.iscomplexobj(y) and y.ndim == 1: + y = N.array([y.real, y.imag]) + else: + raise NotImplementedError,"need complex input for now" + dat = odr.Data(x, y, 1.0 / y**2) + mod = odr.Model(f.fitfcn, extra_args=fcns) + fit = odr.ODR(dat, mod, p0, ifixx=N.zeros(x.ndim), ifixb=fixed, maxit=5000) + fit.run() + return fit.output.beta # should return fit.output + + +