In [ ]:
#%load_ext autoreload
#%autoreload 2
In [ ]:
import os
import bokeh
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from cellpy import cellreader
from cellpy.utils import ica
import holoviews as hv
In [ ]:
import ipywidgets as widgets
from ipywidgets import interact, interact_manual
In [ ]:
%matplotlib inline
hv.extension('bokeh')
In [ ]:
my_data = cellreader.CellpyData()
filename = "../testdata/hdf5/20160805_test001_45_cc.h5"
assert os.path.isfile(filename)
my_data.load(filename)
my_data.set_mass(0.1)
In [ ]:
from lmfit.models import GaussianModel, PseudoVoigtModel, ExponentialGaussianModel, SkewedGaussianModel, LorentzianModel, SkewedVoigtModel, ConstantModel
from lmfit import CompositeModel
The natural way (and my impression is that this is how other groups also do it) of conducting an "in-depth" ica study on a LiB cell would be to measure ica of re-buildt half-cells of both the cathode and the anode of the full cells, fit the peaks of the half cells, and then a use convolution of these fits to fit the actual full cell.
Should include references here...
In [ ]:
class PeakEnsamble:
"""A PeakEnsamble consists of a scale and a set of peaks.
The PeakEnsamble can be fitted with all the internal parameters fixed while only the scale parameter is
varied (jitter=False), or the scale parameter is fixed while the internal parameters (individual peak heights etc)
varied.
Example:
class SiliconPeak(PeakEnsamble):
def __init__():
super().__init__()
self.name = name # Set a prefix name
self.prefixes = [self.name + "Scale", self.name + "01"]]
self.peak_types = [ConstantModel, SkewedGaussianModel]
self.a_new_variable = 12.0
self._read_peak_definitions()
self._init_peaks()
def _read_peak_definitions(self):
self._peak_definitions = {
....
def _init_peaks(self):
self._peaks = self._create_ensamble()
self._set_hints()
self._set_custom_hints()
def _set_custom_hints(self):
....
Attributes:
shift (float): A common shift for all peaks in the ensamble. Should be able to fit
this if jitter is False. TODO: chekc this up.
name (str): Identification label that will be put in front of all peak parameter names.
fixed (bool):
jitter (bool): Allow for individual fitting of the peaks in the ensamble (defaults to True).
max_point (float): The max point of the principal peak.
scale (float): An overall scaling parameter (fitted if jitter=False).
sigma_p1 (float): Sigma value for the principal peak (usually the first peak). When creating the peak
ensamble, parameters for the principal peak is set based on absolute values, while the other
peak parameters are set based on relative values to the principal peak.
prefixes (list): Must be set in the subclass.
"""
def __init__(self, fixed=False, name=None, max_point=1.0,
shift=0.0, sigma_p1=0.01, scale=1.0, jitter=True):
self._peaks = None
self.shift = shift
self.name = name
self.fixed = fixed
self.max_point = max_point
self.jitter = jitter
self.scale = scale
self.sigma_p1 = sigma_p1
self.peak_info = dict()
self._peak_definitions = None
self.peak_var_names = None
self._params = None
def update_peak_parameter(self, name, **kwargs):
self._peaks.set_param_hint(name, **kwargs)
self._make_params()
def reset_peaks(self):
self._init_peaks()
self._custom_init_peaks()
self._make_params()
def init(self):
self._read_peak_definitions()
self._init_peaks()
@property
def peaks(self):
"""lmfit.CompositeModel"""
return self._peaks
@property
def widgets(self):
"""ipywidgets for controlling peak variables"""
raise NotImplementedError
@property
def params(self):
"""lmfit.Parameters (OrderedDict)"""
if self._params is None:
self._make_params()
return self._params
def _make_params(self):
self._params = self._peaks.make_params()
def _read_peak_definitions(self):
raise NotImplementedError("This method must be implemented when sub-classing")
@property
def peak_definitions(self):
return self._peak_definitions
def _init_peaks(self):
self._peaks = self._create_ensamble()
self._set_hints()
def _custom_init_peaks(self):
pass
def _create_ensamble(self):
try:
self.peak_info[self.prefixes[0]] = self.peak_types[0](prefix=self.prefixes[0])
self.peak_info[self.prefixes[1]] = self.peak_types[1](prefix=self.prefixes[1])
except AttributeError:
print("you are missing peak info")
return
p = self.peak_info[self.prefixes[1]]
for prfx, ptype in zip(self.prefixes[2:], self.peak_types[2:]):
self.peak_info[prfx] = ptype(prefix=prfx)
p += self.peak_info[prfx]
p *= self.peak_info[self.prefixes[0]]
return p
def _set_hints(self):
jitter = self.jitter
if self.jitter:
vary = True
vary_scale = False
else:
vary = False
vary_scale = True
scale = self.scale
value_dict = dict()
peak_definitions = self.peak_definitions
prefix_scale = self.prefixes[0]
prefix_peak_1 = self.prefixes[1]
# iterate through all the peaks (not the scale) and collect variables in the value_dict
for var_stub in peak_definitions:
dd = peak_definitions[var_stub]
val_1, ((frac_min, shift_min), (frac_max, shift_max)) = dd[0:2]
v_dict = dict()
v_dict[prefix_peak_1] = [val_1, frac_min * (val_1+shift_min), frac_max * (val_1+shift_max)]
for prfx, (fact, step) in zip(self.prefixes[2:], dd[2:]):
v_dict[prfx] = [fact * (x + step) for x in v_dict[prefix_peak_1]]
value_dict[var_stub] = v_dict
# set parameter hints based on the value_dict
for key1 in value_dict:
for key2 in value_dict[key1]:
_vary = vary
_v = value_dict[key1][key2]
k = "".join((key2, key1))
self._peaks.set_param_hint(k, value=_v[0], min=_v[1], max=_v[2], vary=_vary)
# set parameter hints for scale
self._peaks.set_param_hint("".join((prefix_scale, "c")), value=scale, min=0.1*scale, max=10*scale, vary=vary_scale)
def _fix_full(self, prefix):
"""fixes all variables (but only for this ensamble)"""
for k in self._params:
if k.startswith(prefix):
self._params[k].vary = False
def fit(self, *args, **kwargs):
res = self.peaks.fit(*args, **kwargs)
return res
class Silicon(PeakEnsamble):
"""Peak ensamble for silicon.
This class is a sub-class of PeakEnsamble. Some new attributes are defined
(in addtion to the inhereted attributes).
Attributes:
prefixes (list): A list of peak names used as the prefix when creating the peaks. The firs prefix
should always be for the scale parameter. It is recommended not to play with this attribute.
This attribute is required when subclassing PeakEnsamble
peak_types (list of lmfit peak models): The length of this list must be the same as the length of the
prefixes. It should start with a ConstantModel. This attribute is required when subclassing
PeakEnsamble.
crystalline (bool): Set to true if the Li3.75Si phase exists.
"""
def __init__(self, scale=1.0, crystalline=False, name="Si", max_point=1000, jitter=True,
crystalline_hysteresis=0.0, compress=1.0, expand=1.0,
**kwargs):
"""
Parameters:
scale (float): overall scaling of the peak ensamble
crystalline (bool): set to True if the crystalline peak should be included
name (str): pre-name that will start all parameters
max_point (float): max point
jitter (bool): allow for individual changes between the peaks if True, fix
all individual inter-peak prms if False.
crystalline_hysteresis (float): additional hysteresis for crystalline peak
"""
super().__init__(sigma_p1=0.1, jitter=jitter, scale=scale, max_point=max_point, **kwargs)
self.name = name
self.prefixes = [self.name + x for x in ["Scale", "01", "02", "03"]] # Always start with scale
self.peak_types = [ConstantModel, SkewedGaussianModel, PseudoVoigtModel, PseudoVoigtModel]
self._crystalline = crystalline
self._crystalline_hysteresis = crystalline_hysteresis
self._compress = compress
self._expand = expand
self.init()
def _read_peak_definitions(self):
self._peak_definitions = {
"center": [
0.25 + self.shift, # value
((1.0, -0.1), (1.0, 0.1)), # (frac-min, shift-min), (frac-max, shift-max)
(1.0, 0.21), # (value-frac, value-shift) between peak 1 and peak 2
(1.0, 0.20 + self._crystalline_hysteresis) # (value-frac, value-shift) between peak 1 and peak 3
],
"sigma": [
self.sigma_p1, # value
((0.5*self._expand, 0.0), (2.0*self._compress, 0.0)), # (frac-min, shift-min), (frac-max, shift-max)
(1.0, 0.0), # (value-frac, value-shift) between peak 1 and peak 2
(0.3, 0.0)
],
"amplitude": [
self.sigma_p1 * self.max_point / 0.4, # value
((0.001, 0.0), (100.0, 0.0)), # (frac-min, shift-min), (frac-max, shift-max)
(1.0, 0.0),
(1.0, 0.0)
],
"gamma": [
1.0, # value
((0.001, 0.0), (2.0, 0.0)), # (frac min, shift min), (frac max, shift max)
(1.0, 0.0),
(1.0, 0.0)
]
}
def _custom_init_peaks(self):
self._set_custom_hints()
def _set_custom_hints(self):
if not self._crystalline:
self._unset_crystalline()
def _unset_crystalline(self):
prefix_p3 = self.prefixes[3]
k = "".join([prefix_p3, "amplitude"])
self._peaks.set_param_hint(k, value=0.00001, min=0.000001, vary=False)
k = "".join([prefix_p3, "fraction"])
self._peaks.set_param_hint(k, value=0.00001, min=0.000001, vary=False)
for n in ["center", "sigma"]:
k = "".join([prefix_p3, n])
self._peaks.set_param_hint(k, vary=False)
def _set_crystalline(self):
prefix_p3 = self.prefixes[3]
a = self._peak_definitions["amplitude"][0]
b = self._peak_definitions["amplitude"][3][0]
k = "".join([prefix_p3, "amplitude"])
self._peaks.set_param_hint(k, value=a*b, min=0.000001, vary=True)
for n in ["center", "sigma"]:
k = "".join([prefix_p3, n])
self._peaks.set_param_hint(k, vary=True)
@property
def crystalline(self):
return self._crystalline
@crystalline.setter
def crystalline(self, value):
if not value and self._crystalline:
self._unset_crystalline()
if not self._crystalline and value:
self._set_crystalline()
self._crystalline = value
@property
def widgets(self):
print("overrides PeakEnsamble.widgets property")
print("because it is easier to develop this here and then copy it back to the subclass")
# Need a widget set for each parameter
# Should consist of
# value with max and min
# refine checkbox
# Example
# Name: Si_peak01_amplitude
# [min] -----o--- [max] [value]
# Fixed: x
# The ipywidgets should be coupled to the lmfit.prms somehow
# Need a combined set with all the sub-widgets where the individual peak widgets are
# greyed out if jitter is not selected
class Graphite(PeakEnsamble):
def __init__(self, scale=1.0, name="G", jitter=False, **kwargs):
super().__init__(max_point=10000.0, jitter=jitter, **kwargs)
self.name = name
self.sigma_p1 = 0.01
self.vary = False
self.vary_scale = True
self.prefixes = [self.name + x for x in ["Scale", "01"]] # Always start with scale
self.peak_types = [ConstantModel, LorentzianModel]
self.init()
def _read_peak_definitions(self):
self._peak_definitions = {
"center": [
0.16 + self.shift, # value
((1.0, -0.05), (1.0, 0.05)), # (frac-min, shift-min), (frac-max, shift-max)
# (1.0, 0.21), # (value-frac, value-shift) between peak 1 and peak 2
# (1.0, 0.20) # (value-frac, value-shift) between peak 1 and peak 3
],
"sigma": [
self.sigma_p1, # value
((0.4, 0.0), (5.0, 0.0)), # (frac-min, shift-min), (frac-max, shift-max)
# (1.0, 0.0),
# (0.3, 0.0)
],
"amplitude": [
self.sigma_p1 * self.max_point / 0.4, # value
((0.2, 0.0), (2.0, 0.0)), # (frac-min, shift-min), (frac-max, shift-max)
# (1.0, 0.0),
# (0.002, 0.0)
],
}
In [ ]:
class CompositeEnsamble(PeakEnsamble):
def __init__(self, *ensambles, **kwargs):
super().__init__(self, **kwargs)
self.ensamble = list(ensambles)
self._peaks = None
self._params = None
self._join()
def _join(self):
if len(self.ensamble) > 0:
peaks_left = self.ensamble[0].peaks
prefixes_left = self.ensamble[0].prefixes
params_left = self.ensamble[0].params
if len(self.ensamble) > 1:
for ens in self.ensamble[1:]:
peaks_left += ens.peaks
prefixes_left += ens.prefixes
params_left += ens.params
self._peaks = peaks_left
self._params = params_left
self.prefixes = prefixes_left
def reset_peaks(self):
for ens in self.ensamble:
ens.reset_peaks()
self._join()
#self.make_params()
def make_params(self):
prms = self._peaks.make_params()
return prms
In [ ]:
def get_widgets(parameters):
print(parameters)
In [ ]:
def fix(prefix):
_pars = p.make_params()
for k in _pars:
if k.startswith(prefix):
p[k].vary = False
In [ ]:
def fitplot(v, dq, res, group_title="fit", invert_dq=False, invert_res=False, table=False, width=500, height=500, size=8):
if invert_dq:
dq = -dq
i = 1
if invert_res:
i = -1
raw = hv.Points((v, dq), label="raw", group=group_title).opts(
width=width, height=height, size=size, alpha=0.3,
xlabel="Voltage",
ylabel="dQ/dv",
)
prt = {
"init": hv.Curve((v, i * res.init_fit), group=group_title).opts(alpha=0.5),
"best": hv.Curve((v, i * res.best_fit), group=group_title),
}
parts = res.eval_components()
for key in parts:
if not key.endswith("Scale"):
prt[key] = hv.Curve((v, i * parts[key]), group=group_title)
layout = raw * hv.NdOverlay(prt)
if table:
x = res.best_values
variables = list(x.keys())
values = list(x.values())
lim_min = [res.params[k].min for k in variables]
lim_max = [res.params[k].max for k in variables]
vary = [res.params[k].vary for k in variables]
# expr = [res.params[k].expr for k in variables]
fit_values = {
"var": variables,
"val": values,
"min": lim_min,
"max": lim_max,
"vary": vary,
# "expr": expr,
}
df_fit_values = pd.DataFrame(fit_values)
labels_fit_values = hv.Table(df_fit_values).opts(width=700, height=height)
layout = layout + labels_fit_values
return layout
In [ ]:
# setting peak values
# -------------------
# method one (seems a bit convoluted)
#p_si = Silicon()
#p_si.params['Si02sigma'].min = 0
#p_si.peaks.set_param_hint('Si02sigma', min=0.0)
# method two (simpler?)
p_si = Silicon()
p_si.peaks.set_param_hint('Si02sigma', min=2.2)
p_si.peaks.make_params()
p_si.params["Si02sigma"]
In [ ]:
# method two implemented in PeakEnsamble
p_si = Silicon()
print(f"hint: {p_si.peaks.param_hints['Si02sigma']} val: {p_si.params['Si02sigma']}")
p_si.update_peak_parameter('Si02sigma', min=0.02, vary=False)
print(f"hint: {p_si.peaks.param_hints['Si02sigma']} val: {p_si.params['Si02sigma']}")
p_si.reset_peaks()
print(f"hint: {p_si.peaks.param_hints['Si02sigma']} val: {p_si.params['Si02sigma']}")
In [ ]:
# method two implemented in CompositeEnsamble
p_t = CompositeEnsamble(Silicon(), Graphite())
print(f"hint: {p_t._peaks.param_hints['Si02sigma']} val: {p_t.params['Si02sigma']}")
p_t.update_peak_parameter('Si02sigma', min=0.02, vary=False)
print(f"hint: {p_t.peaks.param_hints['Si02sigma']} val: {p_t.params['Si02sigma']}")
p_t.reset_peaks()
print(f"hint: {p_t.peaks.param_hints['Si02sigma']} val: {p_t.params['Si02sigma']}")
In [ ]:
In [ ]:
cha, volt = my_data.get_dcap(2)
v, dq = ica.dqdv(volt, cha)
In [ ]:
dpeaks = CompositeEnsamble(Silicon(shift=-0.1), Graphite(shift=-0.03))
res = dpeaks.fit(-dq, x=v)
layout_d = fitplot(v, dq, res, invert_res=True, group_title="discharge", table=True)
layout_d
In [ ]:
cha, volt = my_data.get_dcap(6)
v, dq = ica.dqdv(volt, cha)
dpeaks = CompositeEnsamble(Silicon(shift=-0.1), Graphite(shift=-0.03))
res = dpeaks.fit(-dq, x=v)
layout_d = fitplot(v, dq, res, invert_res=True, group_title="discharge", table=True)
layout_d
In [ ]:
In [ ]:
In [ ]:
def f(n, p, **kwargs):
# could invent a better name for this function, maybe...
cha, volt = my_data.get_ccap(n)
v, dq = ica.dqdv(volt, cha)
res = p.fit(dq, x=v)
layout = fitplot(v, dq, res, **kwargs)
return res, layout
In [ ]:
d = dict()
r = dict()
peaks = CompositeEnsamble(Silicon(crystalline=True), Graphite())
cycle_numbers = my_data.get_cycle_numbers()
for n in cycle_numbers:
if n == 2:
print("turning off crystalline peak")
peaks.ensamble[0].crystalline = False
try:
print(f"Fitting cycle {n}", end=" ")
r[n], d[n] = f(n, peaks, width=300)
print("-> OK")
except:
print("OH NO, Failed!")
In [ ]:
%%opts Curve [width=300]
NdLayout = hv.NdLayout(d, kdims='cycle').cols(4)
NdLayout
In [ ]:
def extract_parameter_values(res):
x = res.best_values
variables = list(x.keys())
values = list(x.values())
lim_min = [res.params[k].min for k in variables]
lim_max = [res.params[k].max for k in variables]
vary = [res.params[k].vary for k in variables]
# expr = [res.params[k].expr for k in variables]
fit_values = {
"name": variables,
"value": values,
"min": lim_min,
"max": lim_max,
"vary": vary,
# "expr": expr,
}
df_fit_values = pd.DataFrame(fit_values)
df_fit_values.index.name = "parno"
return df_fit_values
In [ ]:
def combine_parameter_values(r):
frames = list()
for n in r:
df = extract_parameter_values(r[n])
df["cycle"] = n
frames.append(df)
combined = pd.concat(frames).reset_index()
return combined
In [ ]:
result = combine_parameter_values(r)
In [ ]:
Si02center = result[result.name=="Si02center"]
In [ ]:
curve = hv.Curve((Si02center.cycle, Si02center.value), ("Cycle"), ("Position (V)"))
points = hv.Scatter((Si02center.cycle, Si02center.value), ("Cycle"), ("Position (V)")).opts(size=12)
curve * points
In [ ]:
p_not = CompositeEnsamble(Silicon(), Graphite())
p_is = CompositeEnsamble(Silicon(crystalline=True), Graphite())
In [ ]:
# -getting data
cha, volt = my_data.get_ccap(1)
v, dq = ica.dqdv(volt, cha)
In [ ]:
res_is = p_is.fit(dq, x=v)
layout_cryst = fitplot(v, dq, res_is, group_title="With crystalline peak", table=True)
layout_cryst
In [ ]:
res_not = p_not.fit(dq, x=v)
layout_not = fitplot(v, dq, res_not, group_title="Without crystalline peak", table=True)
layout_not
In [ ]:
p_not.ensamble[0].crystalline = True
res_not_is = p_not.fit(dq, x=v)
layout_not_is = fitplot(v, dq, res_not_is, group_title="With added crystalline peak", table=True)
layout_not_is
In [ ]:
In [ ]:
In [ ]:
# -getting data
cha, volt = my_data.get_ccap(2)
v2, dq2 = ica.dqdv(volt, cha)
In [ ]:
p_is.ensamble[0].crystalline = False
# remark! should find out a way to safely transfer or copy the objects if we would like to keep the initial object unchanged
# (this cuould be valuable when using a notebook to prevent problems with excecution order)
res_is2 = p_is.fit(dq2, x=v2)
layout_cryst2 = fitplot(v2, dq2, res_is2, table=True)
layout_cryst2
In [ ]:
# Method one (use the overloaded sum operator to make several peaks)
p1 = Silicon().peaks + Graphite().peaks
pars1 = p1.make_params()
In [ ]:
p1.param_hints
In [ ]:
# Method two (use CompositeEnsamble, the prefered method)
p2 = CompositeEnsamble(Silicon(), Graphite())
In [ ]:
# Method three (consider creating methods for doing this more smoothly)
p3 = CompositeEnsamble()
p3.ensamble.append(Silicon())
p3.ensamble.append(Graphite())
p3._join()
In [ ]:
# Method one and method tow should create the same result (wrt their peaks attribute)
assert p1.param_hints == p2.peaks.param_hints
assert p1.param_hints == p3.peaks.param_hints
In [ ]:
# FITTING CRYSTALLINE PEAK DOES NOT WORK PROPERLY
#p2.ensamble[0].crystalline = True
#p2.ensamble[0].crystalline
In [ ]:
p2.peaks.param_hints
In [ ]:
# initial fit
# -getting data
cha, volt = my_data.get_ccap(1)
v, dq = ica.dqdv(volt, cha)
# -fitting
res = p2.fit(dq, x=v)
print("> OK <".center(80, "-"))
# -creating a plot
layout = fitplot(v, dq, res)
layout
In [ ]:
In [ ]:
# Information about the fit
# -------------------------
# res.data the measurement data (y-values)
# res.init_fit calculated using the initial peak values (y-values)
# res.best_fit the best fit obtained after fitting (y-values)
# res.params contains fitted parameters object (use e.g. in p2.peaks.eval(params=res.params, x=v))
# res.best_values is a dictionary with the best fitted values (i.e a sub-set of res.params)
# res.fit_report() gives a summary of the fit statistics and values
In [ ]:
def f(n, p, **kwargs):
cha, volt = my_data.get_ccap(n)
v, dq = ica.dqdv(volt, cha)
res = p.fit(dq, x=v)
layout = fitplot(v, dq, res, **kwargs)
return layout
In [ ]:
d = dict()
for n in my_data.get_cycle_numbers():
try:
print(f"Fitting cycle {n}", end=" ")
d[n] = f(n, p2, width=300)
print("-> OK")
except:
print("OH NO, Failed!")
In [ ]:
%%opts Curve [width=300]
NdLayout = hv.NdLayout(d, kdims='cycle').cols(4)
NdLayout
In [ ]: