started adding real part, non-lin complex fitting

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 ')
+        f.write("#%7s"%('T'))
-        parfmt = "%.2f" # T formatting
+        parfmt = "%8.2f" # T formatting
         #             if self.Conductivity != None: pass# always true
-        f.write("%8s %8s %8s " % ("e_s", "sig", "pow_sig"))
+        f.write("%9s%9s%9s " % ("e_s", "sig", "pow_sig"))
-        parfmt += " %.3g %.3g %.2f " # conductivity formatting
+        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)
+            f.write("%9s%9s%9s%9s " % enum_peak)
-            parfmt += " %.3g %.3g %.2f %.2f" # peak formatting
+            print enum_peak
-        f.write("high_lim lower_lim") # TODO: store limits
+            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.imag, 'bo', markersize=3, label="Data")
-        pyplot.loglog(self.data.frequency, self.data.epsilon_fit, 'r-', lw=2, label="Fit")
+        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.set_fit_xlimits(10**log10fmin, 10**log10fmax)
-        self.data.fit_limits = [10**log10fmin, 10**log10fmax,ymin,ymax]
         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)
 

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.fit_limits = [f.min(), f.max(), e_imag.min(), e_imag.max()]
-        self.data_curve.setData(f,e_imag)
+        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):

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
+    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
+
+
+