"""Cole-Cole decomposition in conductivity formulation
"""
import numpy as np
import NDimInv.plot_helper
import NDimInv.model_template as mt
import lib_dd.base_class as base_class
import lib_dd.starting_parameters as starting_parameters
import lib_dd.plot_stats as plot_stats
# import conductivity formulation
import sip_models.cond.cc as cc_cond
plt, mpl = NDimInv.plot_helper.setup()
[docs]class decomposition_conductivity(
plot_stats._plot_stats,
base_class.integrated_parameters,
starting_parameters.starting_parameters,
mt.model_template):
def __init__(self, settings):
self.data_format = 'cre_cim'
self.D_base_dims = None
required_keys = ('Nd', 'tausel', 'frequencies', 'c')
for key in required_keys:
if key not in settings:
raise Exception('required key not found: {0}'.format(key))
self.frequencies = settings['frequencies']
self.set_settings(settings)
self.cc = cc_cond.cc(self.settings['frequencies'])
self.set_settings(settings)
[docs] def set_settings(self, settings):
""" Set the settings and call necessary functions
Parameters
----------
settings: dict
contains settings of the decomposition kernel
"""
self.settings = settings
# extract some variables
self.frequencies = self.settings['frequencies']
self.omega = 2.0 * np.pi * self.frequencies
# set min/max tau values corresponding to data limits
self.tau_data_min = 1 / (2 * np.pi * self.frequencies.max())
self.tau_data_max = 1 / (2 * np.pi * self.frequencies.min())
self.tau, self.s, self.tau_f_values = base_class.determine_tau_range(
settings
)
[docs] def convert_parameters(self, pars):
r""" Convert parameters given as (:math:`\sigma_\infty, m_i`) to the
parameterization used by this class.
"""
pars_converted = np.empty_like(pars)
pars_converted[:] = np.log10(pars[:])
return pars_converted
[docs] def convert_pars_back(self, pars):
r""" Convert parameters given in this parameterization back to the
linear state
Here: From :math:`log_{10}(\sigma_\infty), log_{10}(m_i)`
"""
pars_converted = np.empty_like(pars)
pars_converted[:] = 10 ** pars[:]
return pars_converted
def _get_full_pars(self, pars_dec):
# prepare Cole-Cole parameters
sigmai = 10**pars_dec[0][np.newaxis]
m = 10 ** pars_dec[1:]
tau = self.tau
if m.size != tau.size:
raise Exception('m and tau have different sizes!')
c = np.ones_like(m) * self.settings['c']
pars = np.hstack((sigmai, m, tau, c))
return pars
[docs] def forward(self, pars_dec):
"""Forward response of this model
Parameters
----------
pars_dec: or numpy.ndarray
[log10(sigma_infty), log10(m_i)]
Returns
-------
remim: Nx2 numpy.ndarray
with N the nr of frequencies, and the real and the negative
imaginary parts on the second axis
"""
pars = self._get_full_pars(pars_dec)
response = self.cc.response(pars)
creim = response.cre_cim
return creim
[docs] def Jacobian(self, pars_dec):
"""
Parameters
----------
pars_dec: numpy.ndarray
array containing (log10(sigma_infty), log10(m_i)
Returns
-------
J: (2N) X K numpy.ndarray
containing derivatives.
"""
pars = self._get_full_pars(pars_dec)
# real part
real_J = np.concatenate(
(
self.cc.dre_dlog10sigmai(pars)[:, np.newaxis],
self.cc.dre_dlog10m(pars)
),
axis=1
)
imag_J = -np.concatenate(
(
self.cc.dim_dlog10sigmai(pars)[:, np.newaxis],
self.cc.dim_dlog10m(pars)
),
axis=1
)
J = np.concatenate((real_J, imag_J), axis=0).squeeze()
return J
[docs] def get_data_base_dimensions(self):
"""
Return a dict with a description of the data base dimensions. In this
case we have frequencies and re/im data
"""
if self.D_base_dims is None:
self.D_base_dims = {
0: ['frequency', self.frequencies.size],
1: ['cre_cim', 2]
}
return self.D_base_dims
[docs] def get_data_base_size(self):
"""
Return size of flattened base dimensions
"""
return self.frequencies.size * 2
[docs] def get_model_base_dimensions(self):
"""
Return a dict with a description of the model base dimensions. In this
case we have one dimension: the DD parameters (rho0, mi) where m_i
denotes all chargeability values corresponding to the relaxation times.
"""
M_base_dims = {
0: ['rho0_mi', self.tau.size + 1]
}
return M_base_dims
[docs] def compute_par_stats(self, pars):
r"""For a given parameter set (i.e. a fit result), compute relevant
statistical values such as :math:`m_{tot}`, :math:`m_{tot}^n`,
:math:`\tau_{50}`, :math:`\tau_{mean}`, :math:`\tau_{peak}`
This is the way to compute any secondary results based on the fit
results.
Store in self.stat_pars = dict()
"""
base_class.integrated_parameters.compute_par_stats(self, pars)
# self.stat_pars = {}
# the statistical parameters as computed above relate to the
# resistivity formulation. We must correct some of them and add a few
# parameters.
self.stat_pars['sigma_infty'] = self.stat_pars['rho0'].copy()
def sigma0_linear(pars, tau, s, stat_pars):
r"""Compute :math:`sigma0` using math:`\sigma_\infty` and
:math:`m_{tot}`:
.. math::
\sigma_0 = \sigma_\infty \cdot (1 - m_{tot})
"""
sigma0 = 10 ** stat_pars['sigma_infty'] *\
(1 - 10 ** stat_pars['m_tot'])
return sigma0
def sigma0(pars, tau, s, stat_pars):
return np.log10(sigma0_linear(pars, tau, s, stat_pars))
self.stat_pars['sigma0'] = sigma0(pars, self.tau, self.s,
self.stat_pars)
# rho0 is stored in log10, change sign for 1/rho0
self.stat_pars['rho0'] = self.stat_pars['sigma0'] * -1
def mtotn(pars, tau, s, stat_pars):
r"""
Compute the conductivity mtotn:
.. math::
m_{tot}^n = \frac{m_{tot}}{\sigma_0}
"""
mtotn = stat_pars['m_tot'] - stat_pars['rho0']
return mtotn
self.stat_pars['m_tot_n'] = mtotn(pars, self.tau, self.s,
self.stat_pars)
return self.stat_pars