* moved math functions (fit, hn, etc.) to mathlib.py
* fitresults are stored in a better format
This commit is contained in:
parent
5f49684e01
commit
8fc316f2ff
|
|
@ -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 = []
|
||||
258
QDS.py
258
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)
|
||||
|
|
|
|||
98
QDSToolbar.py
Normal file
98
QDSToolbar.py
Normal file
|
|
@ -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
|
||||
113
mathlib.py
Normal file
113
mathlib.py
Normal file
|
|
@ -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
|
||||
Loading…
Reference in New Issue
Block a user