complex fit and stuff

This commit is contained in:
dominik 2022-03-30 17:27:02 +02:00
parent 4ad82cf5b2
commit e19d32b736
8 changed files with 155 additions and 92 deletions

View File

@ -65,7 +65,7 @@ class MultiModel:
self.bounds = [] self.bounds = []
self._kwargs_right = {} self._kwargs_right = {}
self._kwargs_left = {} self._kwargs_left = {}
self._fun_kwargs = {} self.fun_kwargs = {}
# mapping kwargs to kwargs of underlying functions # mapping kwargs to kwargs of underlying functions
self._ext_int_kw = {} self._ext_int_kw = {}
@ -92,7 +92,7 @@ class MultiModel:
if isinstance(func, MultiModel): if isinstance(func, MultiModel):
strcnt = '' strcnt = ''
kw_dict.update(func.fun_kwargs) kw_dict.update(func.fun_kwargs)
self._fun_kwargs.update({k: v for k, v in kw_dict.items()}) self.fun_kwargs.update({k: v for k, v in kw_dict.items()})
self._ext_int_kw.update({k: k for k in kw_dict.keys()}) self._ext_int_kw.update({k: k for k in kw_dict.keys()})
else: else:
@ -102,7 +102,7 @@ class MultiModel:
for k, v in temp_dic.items(): for k, v in temp_dic.items():
key_ = f'{k}_{idx}' key_ = f'{k}_{idx}'
kw_dict[key_] = v kw_dict[key_] = v
self._fun_kwargs[key_] = v self.fun_kwargs[key_] = v
self._ext_int_kw[key_] = k self._ext_int_kw[key_] = k
strcnt = f'({idx})' strcnt = f'({idx})'
@ -116,13 +116,27 @@ class MultiModel:
self.bounds.extend([(None, None)]*len(func.params)) self.bounds.extend([(None, None)]*len(func.params))
def _left_arguments(self, *args, **kwargs): def _left_arguments(self, *args, **kwargs):
kw_left = {k_int: kwargs[k_ext] for k_ext, k_int in self._ext_int_kw.items() if k_ext in self._kwargs_left} kw_left = {}
for k_ext, k_int in self._ext_int_kw.items():
if k_ext in self._kwargs_left:
if not k_ext.startswith('complex_mode'):
kw_left[k_int] = kwargs[k_ext]
else:
kw_left['complex_mode'] = kwargs['complex_mode']
pl = args[:self._param_left] pl = args[:self._param_left]
return pl, kw_left return pl, kw_left
def _right_arguments(self, *args, **kwargs): def _right_arguments(self, *args, **kwargs):
kw_right = {k_int: kwargs[k_ext] for k_ext, k_int in self._ext_int_kw.items() if k_ext in self._kwargs_right} kw_right = {}
for k_ext, k_int in self._ext_int_kw.items():
if k_ext in self._kwargs_right:
if not k_ext.startswith('complex_mode'):
kw_right[k_int] = kwargs[k_ext]
else:
kw_right['complex_mode'] = kwargs['complex_mode']
pr = args[self._param_left:self._param_len] pr = args[self._param_left:self._param_len]
return pr, kw_right return pr, kw_right
@ -142,10 +156,6 @@ class MultiModel:
def right_func(self, x, *args, **kwargs): def right_func(self, x, *args, **kwargs):
return self._right.func(x, *args, **kwargs) return self._right.func(x, *args, **kwargs)
@property
def fun_kwargs(self):
return self._fun_kwargs
def subs(self, x, *args, **kwargs): def subs(self, x, *args, **kwargs):
""" Iterator over all sub-functions (depth-first and left-to-right) """ """ Iterator over all sub-functions (depth-first and left-to-right) """
pl, kw_left = self._left_arguments(*args, **kwargs) pl, kw_left = self._left_arguments(*args, **kwargs)
@ -159,3 +169,18 @@ class MultiModel:
yield from self._right.subs(x, *pr, **kw_right) yield from self._right.subs(x, *pr, **kw_right)
else: else:
yield self._right.func(x, *pr, **kw_right) yield self._right.func(x, *pr, **kw_right)
def sub_name(self):
if isinstance(self._left, MultiModel):
yield from self._left.sub_name()
elif hasattr(self._left, 'name'):
yield self._left.name
else:
yield self.name + '(lhs)'
if isinstance(self._right, MultiModel):
yield from self._right.sub_name()
elif hasattr(self._right, 'name'):
yield self._right.name
else:
yield self.name + '(rhs)'

View File

@ -61,7 +61,7 @@ class FitRoutine(object):
self.result.pop(idx) self.result.pop(idx)
except ValueError: except ValueError:
raise IndexError('Data {} not found'.format(data)) raise IndexError(f'Data {data} not found')
def set_model(self, func, *args, idx=None, **kwargs): def set_model(self, func, *args, idx=None, **kwargs):
if isinstance(func, Model): if isinstance(func, Model):

View File

@ -1,3 +1,5 @@
from __future__ import annotations
import inspect import inspect
from typing import Sized from typing import Sized
@ -61,11 +63,11 @@ class Model(object):
self._int_func = model.func self._int_func = model.func
if hasattr(model, 'subs'): if hasattr(model, 'subs'):
self._int_iter = model.subs self._int_iter = model.subs
self._iter_name = model.sub_name
self.is_multi = True self.is_multi = True
else: else:
self._int_iter = model.func self._int_iter = model.func
try: try:
self.lb, self.ub = list(zip(*model.bounds)) self.lb, self.ub = list(zip(*model.bounds))
except AttributeError: except AttributeError:
@ -78,16 +80,6 @@ class Model(object):
self.fun_kwargs = {k: v.default for k, v in inspect.signature(model.func).parameters.items() self.fun_kwargs = {k: v.default for k, v in inspect.signature(model.func).parameters.items()
if v.default is not inspect.Parameter.empty} if v.default is not inspect.Parameter.empty}
def set_complex(self, state):
if state not in [None, 'complex', 'real', 'imag']:
raise ValueError('"complex" argument is not None, "complex", "real", "imag"')
self.is_complex = state
if state in ['real', 'imag']:
self._complex_part = state
else:
self._complex_part = False
def set_global_parameter(self, idx, p, var=None, lb=None, ub=None, default_bounds=False): def set_global_parameter(self, idx, p, var=None, lb=None, ub=None, default_bounds=False):
if idx is None: if idx is None:
self.parameter = Parameters() self.parameter = Parameters()
@ -127,26 +119,19 @@ class Model(object):
f = self._int_func(x, *p, *self.fun_args, **kwargs) f = self._int_func(x, *p, *self.fun_args, **kwargs)
if self._complex_part:
if self._complex_part == 'real':
return f.real
else:
return f.imag
return f return f
def sub(self, p, x, **kwargs): def sub(self, p, x, **kwargs):
if not kwargs:
kwargs = self.fun_kwargs
if not self.is_multi: if not self.is_multi:
return [self.func(p, x, **kwargs)] return [self.func(p, x, **kwargs)]
else: else:
if not kwargs:
kwargs = self.fun_kwargs
if self._complex_part:
if self._complex_part == 'real':
return [f.real for f in self._int_iter(x, *p, *self.fun_args, **kwargs)]
else:
return [f.imag for f in self._int_iter(x, *p, *self.fun_args, **kwargs)]
return list(self._int_iter(x, *p, *self.fun_args, **kwargs)) return list(self._int_iter(x, *p, *self.fun_args, **kwargs))
def sub_name(self):
if not self.is_multi:
return [self.name]
else:
return list(self._iter_name())

View File

@ -62,11 +62,11 @@ class FitResultCreator:
part_functions = [] part_functions = []
if model.is_multi: if model.is_multi:
for sub_y in model.sub(p_final, _x, **fun_kwargs): for sub_name, sub_y in zip(model.sub_name(), model.sub(p_final, _x, **fun_kwargs)):
if np.iscomplexobj(sub_y): if np.iscomplexobj(sub_y):
part_functions.append(Signal(_x, sub_y)) part_functions.append(Signal(_x, sub_y, name=sub_name))
else: else:
part_functions.append(Points(_x, sub_y)) part_functions.append(Points(_x, sub_y, name=sub_name))
_y = model.func(p_final, _x, **fun_kwargs) _y = model.func(p_final, _x, **fun_kwargs)
resid = model.func(p_final, x_orig, **fun_kwargs) - y_orig resid = model.func(p_final, x_orig, **fun_kwargs) - y_orig
@ -98,7 +98,7 @@ class FitResultCreator:
FitResult(_x, _y, x_orig, y_orig, parameters, fun_kwargs, resid, FitResult(_x, _y, x_orig, y_orig, parameters, fun_kwargs, resid,
nobs, nvar, model.name, stats, nobs, nvar, model.name, stats,
idx=idx, corr=correlation, pcorr=partial_correlation, idx=idx, corr=correlation, pcorr=partial_correlation,
islog=islog, iscomplex=model.is_complex), islog=islog),
part_functions, part_functions,
) )
@ -139,7 +139,7 @@ class FitResultCreator:
class FitResult(Points): class FitResult(Points):
def __init__(self, x, y, x_data, y_data, params, fun_kwargs, resid, nobs, nvar, name, stats, def __init__(self, x, y, x_data, y_data, params, fun_kwargs, resid, nobs, nvar, name, stats,
idx=None, corr=None, pcorr=None, islog=False, iscomplex=None, idx=None, corr=None, pcorr=None, islog=False,
**kwargs): **kwargs):
self.parameter, name = self._prepare_names(params, name) self.parameter, name = self._prepare_names(params, name)
@ -155,7 +155,7 @@ class FitResult(Points):
self.correlation = corr self.correlation = corr
self.partial_correlation = pcorr self.partial_correlation = pcorr
self.islog = islog self.islog = islog
self.iscomplex = iscomplex self.iscomplex = np.iscomplexobj(self.y)
self.x_data = x_data self.x_data = x_data
self.y_data = y_data self.y_data = y_data
self._model_name = name self._model_name = name

View File

@ -6,6 +6,7 @@ from operator import add
from string import ascii_letters from string import ascii_letters
from typing import Dict, List, Tuple from typing import Dict, List, Tuple
import numpy as np
from pyqtgraph import mkPen from pyqtgraph import mkPen
from .fit_forms import FitTableWidget from .fit_forms import FitTableWidget
@ -284,10 +285,11 @@ class QFitDialog(QtWidgets.QWidget, Ui_FitDialog):
parameter['func'] = func parameter['func'] = func
parameter['order'] = order parameter['order'] = order
parameter['len'] = param_len parameter['len'] = param_len
if self._complex[k] is None: parameter['complex'] = self._complex[k]
parameter['complex'] = self._complex[k] if self._complex[k] is not None:
else: for p_k, p_v in parameter['parameter'].items():
parameter['complex'] = ['complex', 'real', 'imag'][self._complex[k]] p_v[1].update({'complex_mode': self._complex[k]})
parameter['parameter'][p_k] = p_v[0], p_v[1]
func_dict[k] = parameter func_dict[k] = parameter
@ -409,23 +411,29 @@ class QFitDialog(QtWidgets.QWidget, Ui_FitDialog):
color = model['color'] color = model['color']
for p, kwargs in parameters.values(): for p, kwargs in parameters.values():
y = f.func(x, *p, **kwargs) if is_complex is not None:
if is_complex is None: y = f.func(x, *p, complex_mode=is_complex, **kwargs)
if np.iscomplexobj(y):
self.preview_lines.append(PlotItem(x=x, y=y.real, pen=mkPen(width=3)))
self.preview_lines.append(PlotItem(x=x, y=y.imag, pen=mkPen(width=3)))
else:
self.preview_lines.append(PlotItem(x=x, y=y, pen=mkPen(width=3)))
else:
y = f.func(x, *p, **kwargs)
self.preview_lines.append(PlotItem(x=x, y=y, pen=mkPen(width=3))) self.preview_lines.append(PlotItem(x=x, y=y, pen=mkPen(width=3)))
if is_complex in [0, 1]:
self.preview_lines.append(PlotItem(x=x, y=y.real, pen=mkPen(width=3)))
if is_complex in [0, 2]:
self.preview_lines.append(PlotItem(x=x, y=y.imag, pen=mkPen(width=3)))
if isinstance(f, MultiModel): if isinstance(f, MultiModel):
for i, s in enumerate(f.subs(x, *p, **kwargs)): sub_kwargs = kwargs.copy()
if is_complex is not None:
sub_kwargs.update({'complex_mode': is_complex})
for i, s in enumerate(f.subs(x, *p, **sub_kwargs)):
pen_i = mkPen(QtGui.QColor.fromRgbF(*color[i])) pen_i = mkPen(QtGui.QColor.fromRgbF(*color[i]))
if is_complex is None: if np.iscomplexobj(s):
self.preview_lines.append(PlotItem(x=x, y=s, pen=pen_i))
if is_complex in [0, 1]:
self.preview_lines.append(PlotItem(x=x, y=s.real, pen=pen_i)) self.preview_lines.append(PlotItem(x=x, y=s.real, pen=pen_i))
if is_complex in [0, 2]:
self.preview_lines.append(PlotItem(x=x, y=s.imag, pen=pen_i)) self.preview_lines.append(PlotItem(x=x, y=s.imag, pen=pen_i))
else:
self.preview_lines.append(PlotItem(x=x, y=s, pen=pen_i))
return self.preview_lines return self.preview_lines

View File

@ -41,17 +41,31 @@ class QFitResult(QtWidgets.QDialog, Ui_Dialog):
self._opts = [(False, False) for _ in range(len(self._results))] self._opts = [(False, False) for _ in range(len(self._results))]
self.residplot = self.graphicsView.addPlot(row=0, col=0) self.residplot = self.graphicsView.addPlot(row=0, col=0)
self.resid_graph = PlotItem(x=[], y=[], symbol='o', symbolPen=None, symbolBrush=mkBrush(color='r'), pen=None) self.resid_graph = PlotItem(x=[], y=[],
symbol='o', symbolPen=None, symbolBrush=mkBrush(color=(174, 199, 232)),
pen=None)
self.resid_graph_imag = PlotItem(x=[], y=[],
symbol='s', symbolPen=None, symbolBrush=mkBrush(color=(255, 127, 14)),
pen=None)
self.residplot.addItem(self.resid_graph) self.residplot.addItem(self.resid_graph)
self.residplot.addItem(self.resid_graph_imag)
self.residplot.setLabel('left', 'Residual') self.residplot.setLabel('left', 'Residual')
self.fitplot = self.graphicsView.addPlot(row=1, col=0) self.fitplot = self.graphicsView.addPlot(row=1, col=0)
self.data_graph = PlotItem(x=[], y=[], symbol='o', symbolPen=None, symbolBrush=mkBrush(color='r'), pen=None) self.data_graph = PlotItem(x=[], y=[],
symbol='o', symbolPen=None, symbolBrush=mkBrush(color=(174, 199, 232)),
pen=None)
self.data_graph_imag = PlotItem(x=[], y=[],
symbol='s', symbolPen=None, symbolBrush=mkBrush(color=(255, 127, 14)),
pen=None)
self.fitplot.addItem(self.data_graph) self.fitplot.addItem(self.data_graph)
self.fitplot.addItem(self.data_graph_imag)
self.fitplot.setLabel('left', 'Function') self.fitplot.setLabel('left', 'Function')
self.fit_graph = PlotItem(x=[], y=[]) self.fit_graph = PlotItem(x=[], y=[])
self.fit_graph_imag = PlotItem(x=[], y=[])
self.fitplot.addItem(self.fit_graph) self.fitplot.addItem(self.fit_graph)
self.fitplot.addItem(self.fit_graph_imag)
self.cmap = RdBuCMap(vmin=-1, vmax=1) self.cmap = RdBuCMap(vmin=-1, vmax=1)
@ -138,15 +152,20 @@ class QFitResult(QtWidgets.QDialog, Ui_Dialog):
res = self._results[idx] res = self._results[idx]
iscomplex = res.iscomplex iscomplex = res.iscomplex
self.resid_graph.setData(x=res.x_data, y=res.residual) if iscomplex:
if iscomplex == 'complex': self.data_graph.setData(x=res.x_data, y=res.y_data.real)
self.data_graph.setData(x=r_[res.x_data, res.x_data], self.data_graph_imag.setData(x=res.x_data, y=res.y_data.imag)
y=r_[res.y_data.real, res.y_data.imag]) self.fit_graph.setData(x=res.x, y=res.y.real)
self.fit_graph.setData(x=r_[res.x, res.x], self.fit_graph_imag.setData(x=res.x, y=res.y.imag)
y=r_[res.y.real, res.y.imag]) self.resid_graph.setData(x=res.x_data, y=res.residual.real)
self.resid_graph_imag.setData(x=res.x_data, y=res.residual.imag)
else: else:
self.resid_graph.setData(x=res.x_data, y=res.residual)
self.resid_graph_imag.setData(x=[], y=[])
self.data_graph.setData(x=res.x_data, y=res.y_data) self.data_graph.setData(x=res.x_data, y=res.y_data)
self.data_graph_imag.setData(x=[], y=[])
self.fit_graph.setData(x=res.x, y=res.y) self.fit_graph.setData(x=res.x, y=res.y)
self.fit_graph_imag.setData(x=[], y=[])
self.fitplot.setLogMode(x=res.islog) self.fitplot.setLogMode(x=res.islog)
self.residplot.setLogMode(x=res.islog) self.residplot.setLogMode(x=res.islog)
@ -243,7 +262,7 @@ class QFitResult(QtWidgets.QDialog, Ui_Dialog):
self.redoFit.emit(self._results) self.redoFit.emit(self._results)
elif button_type == self.buttonBox.Ok: elif button_type == self.buttonBox.Ok:
graph = None graph = '-1'
if self.parameter_checkbox.isChecked(): if self.parameter_checkbox.isChecked():
if self.graph_checkBox.checkState() == QtCore.Qt.Checked: if self.graph_checkBox.checkState() == QtCore.Qt.Checked:
graph = '' graph = ''

View File

@ -390,16 +390,15 @@ class UpperManagement(QtCore.QObject):
models[model_id] = m models[model_id] = m
m_complex = model_p['complex'] m_complex = model_p['complex']
m.set_complex(m_complex)
for set_id, set_params in model_p['parameter'].items(): for set_id, set_params in model_p['parameter'].items():
data_i = self.data[set_id] data_i = self.data[set_id]
if we == 'Deltay': if we == 'Deltay':
we = data_i.y_err**2 we = data_i.y_err**2
if m_complex is None or m_complex == 'real': if m_complex is None or m_complex == 1:
_y = data_i.y.real _y = data_i.y.real
elif m_complex == 'imag' and np.iscomplexobj(self.data[set_id].y): elif m_complex == 2 and np.iscomplexobj(self.data[set_id].y):
_y = data_i.y.imag _y = data_i.y.imag
else: else:
_y = data_i.y _y = data_i.y
@ -518,14 +517,15 @@ class UpperManagement(QtCore.QObject):
if k in parts and show_fit: if k in parts and show_fit:
for subfunc, col in zip(parts[k], TUColorsC): for subfunc, col in zip(parts[k], TUColorsC):
sub_f_id = self.add(subfunc, color=col, linestyle=LineStyle.Dashed, symbol=SymbolStyle.No)
subfunc.value = data_k.value subfunc.value = data_k.value
subfunc.group = data_k.group subfunc.group = data_k.group
sub_f_id = self.add(subfunc, color=col, linestyle=LineStyle.Dashed, symbol=SymbolStyle.No)
f_id_list.append(sub_f_id) f_id_list.append(sub_f_id)
self.delete_sets(tobedeleted) self.delete_sets(tobedeleted)
if accepted and param_graph is not None: if accepted and (param_graph != '-1'):
self.make_fit_parameter(accepted, graph_id=param_graph) self.make_fit_parameter(accepted, graph_id=param_graph)
self.newData.emit(f_id_list, gid) self.newData.emit(f_id_list, gid)
@ -698,7 +698,7 @@ class UpperManagement(QtCore.QObject):
@QtCore.pyqtSlot() @QtCore.pyqtSlot()
def update_color(self): def update_color(self):
UpperManagement._colors = cycle(Colors) UpperManagement._colors = cycle(TUColors)
for i in self.active: for i in self.active:
self.data[i].color = next(UpperManagement._colors) self.data[i].color = next(UpperManagement._colors)
@ -1103,4 +1103,3 @@ class FitWorker(QtCore.QObject):
res = [e.args] res = [e.args]
success = False success = False
self.finished.emit(res, success) self.finished.emit(res, success)

View File

@ -1,5 +1,3 @@
from typing import List, Optional, Tuple
import numpy as np import numpy as np
from ..distributions import Debye, ColeCole, ColeDavidson, KWW, HavriliakNegami from ..distributions import Debye, ColeCole, ColeDavidson, KWW, HavriliakNegami
@ -16,12 +14,18 @@ class _AbstractBDS:
iscomplex = True iscomplex = True
@classmethod @classmethod
def func(cls, x, *args, **kwargs): def func(cls, x, *args, complex_mode: int = 0, **kwargs):
# args[0] : Delta epsilon # args[0] : Delta epsilon
# args[1:] : every other parameter # args[1:] : every other parameter
chi = args[0] * cls.susceptibility(2*np.pi*x, *args[1:]) chi = args[0] * cls.susceptibility(2*np.pi*x, *args[1:], **kwargs)
if complex_mode == 0:
return chi return chi
elif complex_mode == 1:
return chi.real
elif complex_mode == 2:
return chi.imag
else:
raise ValueError(f'{complex_mode!r} is not 0, 1, 2')
class DebyeBDS(_AbstractBDS): class DebyeBDS(_AbstractBDS):
@ -70,9 +74,16 @@ class EpsInfty:
iscomplex = True iscomplex = True
@staticmethod @staticmethod
def func(x, eps): def func(x, eps, complex_mode: int = 0):
ret_val = np.zeros(x.shape, dtype=complex) if complex_mode == 0:
ret_val += eps ret_val = np.zeros(x.shape, dtype=complex)
ret_val += eps
elif complex_mode == 1:
ret_val = eps * np.ones(x.shape)
elif complex_mode == 2:
ret_val = np.zeros(x.shape)
else:
raise ValueError(f'{complex_mode!r} is not 0, 1, 2')
return ret_val return ret_val
@ -86,8 +97,17 @@ class PowerLawBDS:
iscomplex = True iscomplex = True
@staticmethod @staticmethod
def func(x, a, n): def func(x, a, n, complex_mode: int = 0):
return a / (1j*x)**n if complex_mode == 0:
ret_val = np.exp(1j*n*np.pi/2) * a / x**n
elif complex_mode == 1:
ret_val = np.cos(n*np.pi/2) * a / x**n
elif complex_mode == 2:
ret_val = np.sin(n*np.pi/2) * a / x**n
else:
raise ValueError(f'{complex_mode!r} is not 0, 1, 2')
return ret_val
class DCCondBDS: class DCCondBDS:
@ -99,14 +119,21 @@ class DCCondBDS:
iscomplex = True iscomplex = True
@staticmethod @staticmethod
def func(x, sigma): def func(x, sigma, complex_mode: int = 0):
ret_val = np.zeros(x.shape, dtype=complex) if complex_mode == 0:
ret_val += 1j * sigma / x / epsilon0 ret_val = np.zeros(x.shape, dtype=complex)
ret_val += 1j * sigma / x / epsilon0
elif complex_mode == 1:
ret_val = np.zeros(x.shape)
elif complex_mode == 2:
ret_val = sigma / x / epsilon0
else:
raise ValueError(f'{complex_mode!r} is not 0, 1, 2')
return ret_val return ret_val
class HavriliakNegamiDerivative: class DerivativeHavriliakNegami:
name = 'Derivative HN' name = 'Derivative HN'
type = 'Dielectric Spectroscopy' type = 'Dielectric Spectroscopy'
params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma'] params = [r'\Delta\epsilon', r'\tau', r'\alpha', r'\gamma']
@ -123,7 +150,7 @@ class HavriliakNegamiDerivative:
return eps * np.pi * numer / denom / 2. return eps * np.pi * numer / denom / 2.
class ColeColeDerivative: class DerivativeColeCole:
name = 'Derivative CC' name = 'Derivative CC'
type = 'Dielectric Spectroscopy' type = 'Dielectric Spectroscopy'
params = [r'\Delta\epsilon', r'\tau', r'\alpha'] params = [r'\Delta\epsilon', r'\tau', r'\alpha']
@ -140,7 +167,7 @@ class ColeColeDerivative:
return eps * np.pi * numer / denom / 2. return eps * np.pi * numer / denom / 2.
class ColeDavidsonDerivative: class DerivativeColeDavidson:
name = 'Derivative CD' name = 'Derivative CD'
type = 'Dielectric Spectroscopy' type = 'Dielectric Spectroscopy'
params = [r'\Delta\epsilon', r'\tau', r'\gamma'] params = [r'\Delta\epsilon', r'\tau', r'\gamma']
@ -149,7 +176,7 @@ class ColeDavidsonDerivative:
@staticmethod @staticmethod
def func(x, eps, tau, g): def func(x, eps, tau, g):
omtau = 2*np.pi*x * tau omtau = 2*np.pi*x * tau
numer = g * omtau * np.sin((1+g)*np.sin(omtau)) numer = g * omtau * np.sin((1+g)*np.arctan(omtau))
denom = (1 + omtau**2)**((1+g)/2.) denom = (1 + omtau**2)**((1+g)/2.)
return eps * np.pi * numer / denom / 2. return eps * np.pi * numer / denom / 2.