# Source code for impedance.validation

import numpy as np
from impedance.models.circuits.fitting import rmse
from impedance.models.circuits.elements import circuit_elements, K  # noqa

""" A method for implementing the Lin-KK test for validating linearity 

Parameters
----------
f: np.ndarray
measured frequencies
Z: np.ndarray of complex numbers
measured impedances
c: np.float
cutoff for mu
max_M: int
the maximum number of RC elements
fit_type: str
selects which components of data are fit ('real', 'imag', or
'complex')
option to add a serial capacitance that helps validate data with no
low-frequency intercept

Returns
-------
M: int
number of RC elements used
mu: np.float
under- or over-fitting measure
Z_fit: np.ndarray of complex numbers
impedance of fit at input frequencies
resids_real: np.ndarray
real component of the residuals of the fit at input frequencies
resids_imag: np.ndarray
imaginary component of the residuals of the fit at input frequencies

Notes
-----

The lin-KK method from Schönleber et al.  is a quick test for checking
the
validity of EIS data. The validity of an impedance spectrum is analyzed by
its reproducibility by a Kramers-Kronig (KK) compliant equivalent circuit.
In particular, the model used in the lin-KK test is an ohmic resistor,
:math:R_{Ohm}, and :math:M RC elements.

.. math::

\\hat Z = R_{Ohm} + \\sum_{k=1}^{M} \\frac{R_k}{1 + j \\omega \\tau_k}

The :math:M time constants, :math:\\tau_k, are distributed
logarithmically,

.. math::
\\tau_1 = \\frac{1}{\\omega_{max}} ; \\tau_M = \\frac{1}{\\omega_{min}}
; \\tau_k = 10^{\\log{(\\tau_{min}) + \\frac{k-1}{M-1}\\log{{(
\\frac{\\tau_{max}}{\\tau_{min}}}})}}

and are not fit during the test (only :math:R_{Ohm} and :math:R_{k}
are free parameters).

In order to prevent under- or over-fitting, Schönleber et al. propose using
the ratio of positive resistor mass to negative resistor mass as a metric
for finding the optimal number of RC elements.

.. math::

\\mu = 1 - \\frac{\\sum_{R_k \\ge 0} |R_k|}{\\sum_{R_k < 0} |R_k|}

The argument :code:c defines the cutoff value for :math:\\mu. The
algorithm starts at :code:M = 3 and iterates up to :code:max_M until a
:math:\\mu < c is reached. The default of 0.85 is simply a heuristic
value based off of the experience of Schönleber et al., but a lower value
may give better results.

If the argument :code:c is :code:None, then the automatic determination
of RC elements is turned off and the solution is calculated for
:code:max_M RC elements. This manual mode should be used with caution as
under- and over-fitting should be avoided.

 Schönleber, M. et al. A Method for Improving the Robustness of
linear Kramers-Kronig Validity Tests. Electrochimica Acta 131, 20–27 (2014)
doi: 10.1016/j.electacta.2014.01.034
<https://doi.org/10.1016/j.electacta.2014.01.034>_.

"""

if c is not None:
M = 0
mu = 1
while mu > c and M <= max_M:
M += 1
ts = get_tc_distribution(f, M)

if M % 10 == 0:
print(M, mu, rmse(eval_linKK(elements, ts, f), Z))
else:
M = max_M
ts = get_tc_distribution(f, M)

resids_real = residuals_linKK(elements, ts, Z, f, residuals='real')
resids_imag = residuals_linKK(elements, ts, Z, f, residuals='imag')
return M, mu, Z_fit, resids_real, resids_imag

[docs]def get_tc_distribution(f, M):
""" Returns the distribution of time constants for the linKK method """

t_max = 1/(2 * np.pi * np.min(f))
t_min = 1/(2 * np.pi * np.max(f))
ts = np.zeros(shape=(M,))
ts = t_min
ts[-1] = t_max
if M > 1:
for k in range(2, M):
ts[k-1] = 10**(np.log10(t_min) +
((k-1)/(M-1))*np.log10(t_max/t_min))
return ts

""" Fits the linKK model using linear regression

Parameters
----------
f: np.ndarray
measured frequencies
ts: np.ndarray
logarithmically spaced time constants of RC elements
M: int
the number of RC elements
Z: np.ndarray of complex numbers
measured impedances
fit_type: str
selects which components of data are fit ('real', 'imag', or
'complex')
option to add a serial capacitance that helps validate data with no
low-frequency intercept

Returns
-------
elements: np.ndarray
values of fit :math:R_k in RC elements and series :math:R_0,
L, and optionally C.
mu: np.float
under- or over-fitting measure

Notes
-----
Since we have a system of equations, :math:Ax ~= b, that's linear wrt
:math:R_k, we can fit the model by calculating the pseudo-inverse of A.
:math:Ax is our model fit, :math:\\hat{Z}, and :math:b is the
normalized real or imaginary component of the impedance data,
:math:Re(Z)/|Z| or :math:Im(Z)/|Z|, respectively.

:math:\\hat{Z} = R_0 + \\sum^M_{k=1}(R_k / |Z|(1 + j * w * \\tau_k)).
:math:x is an (M+1) :math:\\times 1 matrix where the first row
contains :math:R_0 and subsequent rows contain :math:R_k values.
A is an N :math:\\times (M+1) matrix, where N is the number of data
points, and M is the number of RC elements.

Examples
--------

Fitting the real part of data, the first column of A contains
values of :math:\\frac{1}{|Z|}, the second column contains
:math:Re(1 / |Z| (1 + j * w * \\tau_1)), the third contains
:math:Re(1 / |Z| (1 + j * w * \\tau_2)) and so on. The :math:R_k values
within the x matrix are found using :code:numpy.linalg.pinv when
fit_type = 'real' or 'imag'. When fit_type = 'complex' the coefficients are
found "manually" using :math:r = ||A'x - b'||^2 + ||A''x - b'||^2
according to Eq 14 of Schonleber .

 Schönleber, M. et al. A Method for Improving the Robustness of
linear Kramers-Kronig Validity Tests. Electrochimica Acta 131, 20–27 (2014)
doi: 10.1016/j.electacta.2014.01.034
<https://doi.org/10.1016/j.electacta.2014.01.034>_.
"""

w = 2 * np.pi * f

# Fitting model has M RC elements plus 1 series resistance and 1 series
# inductance
a_re = np.zeros((f.size, M+2))
a_im = np.zeros((f.size, M+2))

a_re = np.zeros((f.size, M+3))
a_im = np.zeros((f.size, M+3))

# Column for series capacitance. Real part = 0.
a_im[:, -2] = - 1 / (w * np.abs(Z))

# Column for series resistance, R_o in model. Imaginary part = 0.
a_re[:, 0] = 1 / np.abs(Z)

# Column for series inductance to capture inevitable contributions from
# the measurement system. Real part = 0.
a_im[:, -1] = w / np.abs(Z)

# Columns for series RC elements
for i, tau in enumerate(ts):
a_re[:, i+1] = K([1, tau], f).real / np.abs(Z)
a_im[:, i+1] = K([1, tau], f).imag / np.abs(Z)

if fit_type == 'real':
elements = np.linalg.pinv(a_re).dot(Z.real / np.abs(Z))

# After fitting real part, need to use imaginary component of fit to
# find values of series inductance and capacitance
a_im = np.zeros((f.size, 2))
a_im[:, -1] = w / np.abs(Z)
a_im[:, -2] = -1 / (w * np.abs(Z))
elements[-2] = 1e-18  # nullifies series C without dividing by 0

coefs = np.linalg.pinv(a_im).dot((Z.imag - Z_fit_re.imag) / np. abs(Z))

elements[-2:] = coefs
else:
elements[-1] = coefs[-1]
elif fit_type == 'imag':
elements = np.linalg.pinv(a_im).dot(Z.imag / np.abs(Z))

# Calculates real part of impedance from fitting to imaginary parts
# without ohmic resistance, i.e. only the real parts of series RC
# elements.

# Weighting used in Boukamp et al - "A Linear Kronig-Kramers
# Transform for Immittance Data Validation" 1995, J. Electrochem
# Soc. 142 (6)
ws = 1 / (Z.real**2 + Z.imag**2)

# Finds ohmic resistance for imaginary part fit according to Eq 7
# of Boukamp et al.
elements = np.sum(ws * (Z.real - z_re.real)) / np.sum(ws)
elif fit_type == 'complex':
# x = (A*•A)^-1
# y = A*•b
# Pseudoinsverse, A^+ = (A*•A)^-1•A* and R_k values = A^+•b
x = np.linalg.inv(a_re.T.dot(a_re) + a_im.T.dot(a_im))
y = a_re.T.dot(Z.real / np.abs(Z)) + a_im.T.dot(Z.imag / np.abs(Z))
elements = x.dot(y)
else:
raise ValueError("Invalid choice of fit_type, please choose from "
"\'real\', \'imag\', or \'complex\'")

mu = calc_mu(elements[1:-2])
else:
mu = calc_mu(elements[1:-1])

return elements, mu

""" Builds a circuit of RC elements to be used in LinKK """
circuit_string = 's([R({},{}),'.format([elements], f.tolist())

for (Rk, tk) in zip(elements[1:], ts):
circuit_string += f'K({[Rk, tk]},{f.tolist()}),'

circuit_string += 'L({},{}),'.format([elements[-1]], f.tolist())
if elements.size == (ts.size + 3):
circuit_string += 'C({},{}),'.format([1/elements[-2]], f.tolist())

circuit_string = circuit_string.strip(',')
circuit_string += '])'

return eval(circuit_string, circuit_elements)

[docs]def residuals_linKK(elements, ts, Z, f, residuals='real'):
""" Calculates the residual between the data and a LinKK fit """

err = (Z - eval_linKK(elements, ts, f))/np.abs(Z)

if residuals == 'real':
return err.real
elif residuals == 'imag':
return err.imag
elif residuals == 'both':
z1d = np.zeros(Z.size*2, dtype=np.float64)
z1d[0:z1d.size:2] = err.real
z1d[1:z1d.size:2] = err.imag

return z1d

[docs]def calc_mu(Rs):
""" Calculates mu for use in LinKK """

neg_sum = sum(abs(x) for x in Rs if x < 0)
pos_sum = sum(abs(x) for x in Rs if x >= 0)

return 1 - neg_sum/pos_sum