# -*- encoding: utf-8 -*- __author__ = 'markusro' from PyQt4.QtGui import QColor from PyQt4.QtCore import QObject,pyqtSignal,QThread,pyqtSlot import numpy as np from scipy import optimize as opt, odr import libyaff def fit_anneal(x, y, p0, fixed, funcs): raise NotImplementedError 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, funcs): raise NotImplementedError # 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 # Replaced with fit_odr_cmplx # #def fit_odr(x, y, p0, fixed, funcs): # dat = odr.Data(x, y, 1.0 / y**2) # 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 * np.pi * nu Phi = np.arctan((om * tau) ** a * np.sin(np.pi * a / 2.) / (1. + (om * tau) ** a * np.cos(np.pi * a / 2.))) e_loss = delta_eps * (1 + 2 * (om * tau) ** a * np.cos(np.pi * a / 2.) + (om * tau) ** (2. * a) ) ** ( -b / 2.) * np.sin(b * Phi) e_stor = delta_eps * (1 + 2 * (om * tau) ** a * np.cos(np.pi * a / 2.) + (om * tau) ** (2. * a) ) ** ( -b / 2.) * np.cos(b * Phi) return e_loss # 2* oder nicht? def id_to_color(id): colors = [ QColor(54,22,115), QColor(160,16,36), QColor(45,142,15), QColor(168,149,17), QColor(36,10,85), QColor(118,8,23), QColor(31,105,7), QColor(124,109,8), ] return colors[id % len(colors)] def mini_func(p, x, y): res = y - multi_hn(p, x) # apply weights res /= 1 / y return np.sqrt(np.dot(res, res)) def multi_hn(p, nu): conductivity = p[1] cond_beta = p[2] om = 2 * np.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 = np.arctan((om * tau) ** a * np.sin(np.pi * a / 2.) / (1. + (om * tau) ** a * np.cos(np.pi * a / 2.))) e_loss += 2 * delta_eps * (1 + 2 * (om * tau) ** a * np.cos(np.pi * a / 2.) + (om * tau) ** (2. * a) ) ** ( -b / 2.) * np.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 tau_peak(f, a, b): tau = (np.sin(np.pi * a / 2. / (b + 1)) / np.sin(np.pi * a * b / 2. / (b + 1))) ** (1 / a) tau /= 2 * np.pi * f return tau ### define funcs here class Functions(QObject): def __init__(self): super(Functions,self).__init__() self.list = { # provides functions: "id_string":(function, number_of_parameters) "hn":(self.hn_cmplx, 4), "conductivity":(self.cond_cmplx, 3), "power":(self.power_cmplx, 2), "static":(self.static_cmplx, 1), "yaff":(self.yaff, 8) } self.YAFF = libyaff.Yaff() def hn_cmplx(self, p, x): om = 2*np.pi*x #hn = om*1j eps,t,a,b = p hn = eps/(1+(1j*om*t)**a)**b cplx = np.array([hn.real, -hn.imag]) return cplx def cond_cmplx(self, p, x): om = 2*np.pi*x sgma, isgma, n = p cond = sgma/(om**n) + isgma/(1j*om**n) # Jonscher (Universal Dielectric Response: e",e' prop sigma/omega**n cplx = np.array([cond.real, -cond.imag]) return cplx def power_cmplx(self, p, x): om = 2*np.pi*x sgma,n = p power = sgma/(om*1j)**n cplx = np.array([power.real, -power.imag]) return cplx def static_cmplx(self, p, x): eps_inf = p[0] static = np.ones( (2,x.size) )*eps_inf static[1,:] *= 0 # set imag part zero #cplx = N.array([static.real, static.imag]) return static def yaff(self,p,x): ya = self.YAFF.yaff(p[:8],x) cplx = np.array([ya.imag, ya.real]) return cplx def get(self,name): return self.list[name] def get_function(self,name): return self.list[name][0] class FitFunctionCreator(QObject): new_data = pyqtSignal(np.ndarray, np.ndarray) def __init__(self): super(FitFunctionCreator,self).__init__() self.data = None self.functions = Functions() def fitfcn(self, p0, x, *funcs): if x.ndim == 2: self.data = np.zeros( x.shape ) else: self.data = np.zeros( (2,x.size) ) ndx = 0 for fn in funcs: # loop over functions and add the results f, num_p = fn.function, fn.param_number p = p0[ndx:ndx + num_p] if x.ndim == 2: x = x[0] result = f(p, x) self.data += result # fit functions take only 1-dim x ndx += num_p self.new_data.emit(x, self.data) return self.data def _fitfcn(self, p0, x, *funcs): if x.ndim == 2: self.data = np.zeros( x.shape ) else: self.data = np.zeros( (2,x.size) ) ndx = 0 for fn in funcs: # loop over functions and add the results f,num_p = self.functions.get(fn) p = p0[ndx:ndx + num_p] if x.ndim == 2: x = x[0] result = f(p, x) self.data += result # fit functions take only 1-dim x ndx += num_p self.new_data.emit(x, self.data) 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 np.iscomplexobj(y) and y.ndim == 1: weights = 1/np.abs(y)**2 we = np.resize(weights, (2, weights.size)) #we = 1/N.array([y.real**2, y.imag**2]) y = np.array([y.real, y.imag]) else: raise NotImplementedError, "need complex input for now" dat = odr.Data(x, y, we=we) mod = odr.Model(f.fitfcn, extra_args=fcns) fit = odr.ODR(dat, mod, p0, ifixx=(0,), ifixb=fixed, maxit=8000) fit.run() #print fit.output.pprint() return fit.output class FitRoutine(QObject): finished_fit = pyqtSignal() data_ready = pyqtSignal(np.ndarray, np.ndarray) def __init__(self): super(FitRoutine,self).__init__() self.f = FitFunctionCreator() self.f.new_data.connect(self.data_ready.emit) self._fitter = self.fit_odr_cmplx self._odr_fit = None @property def fitter(self): return self._fitter @fitter.setter def fitter(self,f): self._fitter = f def fit_odr_cmplx(self, x, y, p0, fixed, fcns): #if x.ndim < 2: # x = N.resize(x, (2,x.size)) if np.iscomplexobj(y) and y.ndim == 1: weights = 1/np.abs(y)**2 we = np.resize(weights, (2, weights.size)) #we = 1/N.array([y.real**2, y.imag**2]) y = np.array([y.real, y.imag]) else: raise NotImplementedError, "need complex input for now" dat = odr.Data(x, y, we=we) mod = odr.Model(self.f.fitfcn, extra_args=fcns) self._odr_fit = odr.ODR(dat, mod, p0, ifixx=(0,), ifixb=fixed, maxit=10) @pyqtSlot() def fit(self): #print "TID in FitRoutine", QThread.thread() self._odr_fit.run() self.finished_fit.emit() def result(self): return self._odr_fit.output class FunctionRegister: def __init__(self): self.registry = {} def register_function(self, obj): #print "FR: Registering:",obj id_string = obj.id_string if self.registry.has_key(obj): raise AssertionError,"The object is already registered! This should NOT happen" self.registry[obj]=id_string #print "FR: ",self.registry def unregister_function(self, obj): #print "FR: UnRegistering:",obj if self.registry.has_key(obj): self.registry.pop(obj) else: raise AssertionError,"The object is not in the registry! This should NOT happen" #print "FR: ",self.registry def get_registered_functions(self): return self.registry