<jupyter_start><jupyter_text>Machine learning pipelines using battery dataAuthors: Dr. <NAME>, <NAME>, Dr. <NAME> (NREL)This notebook demonstrates how to use machine-learning to develop predictive models for battery health diagnosis, using electrochemical impedance spectroscopy (EIS) to predicy discharge capacity. This challenge is used as an example to demonstrate fundamentals of machine-learning, such as data standardization, train/test splits, feature engineering, and hyperparameter optimization. This tutorial also serves as a summary of how to use the machine-learning library sklearn, especially the functionality of Pipelines, which are extremely useful tools for automating the data processing and model optimization process in a way that is repeatable and understandable.This tutorial uses open-source data from Zhang et al, *Nature Communications* (2020) 11:1706 ([paper](https://www.nature.com/articles/s41467-020-15235-7.pdf), [github](https://github.com/YunweiZhang/ML-identify-battery-degradation)).  Setup<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from eisplot import plot_eis
from plotyy import plotyy
# Load data
eis_data = pd.read_csv("data/data_eis_zhang2020.csv")
freq = np.loadtxt("data/freq_eis_zhang2020.csv")<jupyter_output><empty_output><jupyter_text>Processing impedance dataInvestigate the dataframe.<jupyter_code>eis_data
eis_data.describe()<jupyter_output><empty_output><jupyter_text>There are several variables in this data frame:- seriesIdx: index denoting which cell (1-8) was the data is from- cycle: cycle number count- Q: capacity of the cell in mAh- q: relative capacity of each cell (capacity divided by initial capacity, calculated for each cell)- ZReal_xHz: real component of the impedance at x Hz frequency- ZImag_xHz: imaginary component of the impedance at x Hz frequency- Z_xHz: magnitude of the impedance at x Hz frequency- ZPhase_xHz: phase angle of the impedance at x Hz frequencySo, for training a model, the target variable is Q, while the features are all the components of the measured impedance at each frequency. The frequency vector for these measurements can be created using `np.logspace` or loaded from file.<jupyter_code># Frequency vector
freq.shape, np.max(freq), np.min(freq)
y = eis_data['Q']
is_Z_col = eis_data.columns.str.contains('Z')
x = eis_data[eis_data.columns[is_Z_col]]
y.shape, x.shape<jupyter_output><empty_output><jupyter_text>The data is now extracted, and the shapes look correct (same number of rows in x and y, 240 columns in x corresponds to 4 components of the impedance for 60 frequencies).To plot the impedance data, we need to reform these features into arrays of complex numbers.<jupyter_code>is_ZReal_col = x.columns.str.contains('ZReal')
is_ZImag_col = x.columns.str.contains('ZImag')
eis = x.apply(lambda d: d[is_ZReal_col].values + 1j*d[is_ZImag_col].values, axis=1)
eis[0].shape<jupyter_output><empty_output><jupyter_text>Visualizing the dataVisualizing the raw data can help us understand the complexity of the problem. We should plot both the capacity (our target variable) as well as all of our impedance data (our features). The task of plotting both features and target variables can be very challenging, but is very important for getting an intuitive understanding of your data set. Visualizing capacityPlot the capacity of each cell versus cycle count.<jupyter_code>series = np.unique(eis_data.seriesIdx)
# Customize the colors to match with scatter plots later in the notebook
cmap = plt.get_cmap('tab10', series.size)
for this_series in series:
    mask_series = eis_data.seriesIdx == this_series
    plt.plot(eis_data.cycle[mask_series], eis_data.Q[mask_series], '-',
        color=cmap((this_series-1)/series.size), linewidth=2, label='Series %d' % (this_series))
plt.ylim([0, 40])
plt.xlim([0, 300])
plt.xlabel('Cycle count')
plt.ylabel('Discharge capacity (mAh)')

plt.legend()<jupyter_output><empty_output><jupyter_text>There's wide variation in the magnitude of the initial drop, and then all cells stabilize for between 20 cycles (before testing is cutoff, series 4 and 8) to about 200 cycles (series 1 and 5). Some cells display an extremely sudden drop of capacity at 100 cycles (series 7), while others have just a slight increase in the degradation rate at 100 cycles (series 1 and 6), and other cells have a mix of behaivors (series 1: change of slope at 100 cycles, sudden drop around 200 cycles).  Visualizing impedance dataPlotting just one curve gives us an idea of the typical 'shape' of an impedance measurement, as well as the relationships between the components of the impedance. |Z| and Z_Real are both 'magnitude' features, which are always increasing as frequency decreases, while Z_Imaginary and phase of Z are both 'phase' features, which vary non-monotonically versus frequency.<jupyter_code>plot_eis(freq, eis[0])<jupyter_output><empty_output><jupyter_text>Plotting all of the data for one of the cells can help us understand how capacity and impedance vary throughout the lifetime of a cell.<jupyter_code>this_series = 1
mask_series = eis_data.seriesIdx == this_series

plt.scatter(eis_data.cycle[mask_series], eis_data.Q[mask_series],
            c=eis_data.cycle[mask_series], cmap='viridis')
plt.xlabel('Cycle count')
plt.ylabel('Discharge capacity (mAh)')
plt.ylim([0, 40])
plt.xlim([0, 300])
plt.title('Series %d' % (this_series))

plot_eis(freq, eis[mask_series], cmap='viridis', title='Series %d' % (this_series))<jupyter_output><empty_output><jupyter_text>Impedance seems to be increasing throughout the lifetime of the cell. Also, the locations of peaks and valleys in the phase of Z and Z_Imaginary seem to vary. Let's plot the data from one more cell.<jupyter_code>this_series = 2
mask_series = eis_data.seriesIdx == this_series

plt.scatter(eis_data.cycle[mask_series], eis_data.Q[mask_series],
            c=eis_data.cycle[mask_series], cmap='viridis')
plt.xlabel('Cycle count')
plt.ylabel('Discharge capacity (mAh)')
plt.ylim([0, 40])
plt.xlim([0, 300])
plt.title('Series %d' % (this_series))

plot_eis(freq, eis[mask_series], cmap='viridis', title='Series %d' % (this_series))<jupyter_output><empty_output><jupyter_text>This battery capacity rapidly fades at beginning-of-life before stabilizing at around 26 mAh. Throughout this the impedance decreases slightly.The wide variance between trends implies that it should be challenging to predict capacity from impedance.  Training a model We can train a simple linear model predicting capacity with very few lines of code using sklearn. Best practice is to normalize the input features before training.<jupyter_code># Import regression model class
from sklearn.linear_model import LinearRegression
# Import data normalization tool
from sklearn.preprocessing import StandardScaler
# Import standard statistical metrics for evaluating model performance
from sklearn.metrics import mean_absolute_error, r2_score
# Normalize the training data
scaler = StandardScaler()
x_ = scaler.fit_transform(x)
# Instantiate the linear model class object
regressor = LinearRegression()
# Fit the linear model
regressor.fit(x_, y)
# Make the predictions using the fit model
y_pred = regressor.predict(x_)
print("MAE: %0.2f, R^2: %0.2f" % (mean_absolute_error(y, y_pred), r2_score(y, y_pred)))
plotyy(y, y_pred, color=eis_data.seriesIdx)<jupyter_output><empty_output><jupyter_text>While this looks impressive, it's seems suspiciously good. Should we really be able to predict the capacity of a battery with such high accuracy after making so little effort?One common method to test our model performance is to hold out data from training to test a model on, i.e., the test set. A common way to test a model is to randomly hold out some fraction of the data.<jupyter_code>from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(x_, y,
    test_size=0.3, random_state=42)

x_train.shape, x_test.shape
# Train
regressor.fit(x_train, y_train)
# Test
y_test_pred = regressor.predict(x_test)
print("MAE: %0.2f, R^2: %0.2f" % (mean_absolute_error(y_test, y_test_pred), r2_score(y_test, y_test_pred)))
# Note we can access the seriesIdx using the indices of the test dataframe, as sklearn plays very nicely with pandas
plotyy(y_test, y_test_pred, color=eis_data.seriesIdx[y_test.index])<jupyter_output>MAE: 0.46, R^2: 0.99<jupyter_text>This is only a tiny bit worse. Our simple data exploration seemed to suggest that this should be challenging task.Here, it turns out that our test set is actually not very difficult - we are simply picking points at random, regardless of what cell they are from. In the 'real world', we need to predict the capacity of cells that we have no capacity data for.  Predicting capacity from unseen cellsWhen predicting from new data, we have the additional complication of having to 'train' any data transformations, and then replicated them on the test set. For instance, we should normalize our data only using information from the training set, and then transform the test set data accordingly. This practice prevents 'data leakage', e.g., spurious correlations - we don't want **any** information from the test data to affect our model during training; if we normalize the entire data set at once, we are including statistical information about the test data into the normalization statistics.<jupyter_code># Use the last cell for testing
is_train = eis_data.seriesIdx != eis_data.seriesIdx.iloc[-1]
x_train, x_test, y_train, y_test = x[is_train], x[~is_train], y[is_train], y[~is_train]
# Standardize the training data
scaler = StandardScaler()
x_train_ = scaler.fit_transform(x_train)
# Train
regressor.fit(x_train_, y_train)
# Test
x_test_ = scaler.transform(x_test)
y_test_pred = regressor.predict(x_test_)
title = "Testing on an unseen cell \n" r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y_test, y_test_pred), r2_score(y_test, y_test_pred))
plotyy(y_test, y_test_pred, title=title)<jupyter_output><empty_output><jupyter_text>We could also plot the actual and predicted versus cycle count:<jupyter_code>plt.plot(eis_data.cycle[y_test.index], y_test, 'ok')
plt.plot(eis_data.cycle[y_test.index], y_test_pred, 'dr')
plt.xlabel('Cycle count')
plt.ylabel('Discharge capacity (mAh)')
plt.legend(('Actual', 'Predicted'))<jupyter_output><empty_output><jupyter_text>We should also always be comparing results with some kind of baseline, to contextualize the error statistics. For this relatively straightforward problem, it's not hard that hard to contextualize the results, but for more complex problems, it can be difficult to understand the error statistics. It's also simply best practice to make sure your models are outperforming extremely simple approaches. There are several 'dummy' model types available in sklearn to compare against.<jupyter_code>from sklearn.dummy import DummyRegressor

# Use the last cell for testing
is_train = eis_data.seriesIdx != eis_data.seriesIdx.iloc[-1]
x_train, x_test, y_train, y_test = x[is_train], x[~is_train], y[is_train], y[~is_train]
# Standardize the training data
scaler = StandardScaler()
x_train_ = scaler.fit_transform(x_train)
# Train
dummy_regr = DummyRegressor(strategy='mean') # this is the default setting
dummy_regr.fit(x_train_, y_train)
# Test
x_test_ = scaler.transform(x_test)
y_test_pred = dummy_regr.predict(x_test_)
title = "Dummy model test \n" r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y_test, y_test_pred), r2_score(y_test, y_test_pred))
plotyy(y_test, y_test_pred, title=title)
# Capacity vs. cycle plot
fig, ax = plt.subplots(1,1)
plt.plot(eis_data.cycle[y_test.index], y_test, 'ok')
plt.plot(eis_data.cycle[y_test.index], y_test_pred, 'dr')
plt.xlabel('Cycle count')
plt.ylabel('Discharge capacity (mAh)')
plt.legend(('Actual', 'Predicted'))<jupyter_output><empty_output><jupyter_text>It's pretty clear that this model is no longer working well - the dummy model is actually performing better (MAE of 1.3 mAh, vs. 3 mAh for the linear model). But our test set is quite limited. We want to know how well this model performs on **any** unseen cell, not just **one** unseen cell. For this, we can use cross-validation, where we repeat this process, but holding out each cell in turn for testing and training on the other cells. But recall that we need to 'train' the data standardization for each training set, as well. It's not too hard to imagine how this process would work, but it is quite a few lines of code. Thankfully, sklearn automates this process using 'Pipelines'. Model pipelinesWe can set up a simple example pipeline using the easy train/test split from earlier (randomized splitting by data point). Below is one way to instantiate a pipeline. You can also use the `make_pipeline()` method.<jupyter_code>from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_predict
from sklearn.model_selection import KFold

# Instatiate the pipeline object
pipe = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('regressor', LinearRegression())
    ]
)
# Run cross-validation
# Default is 5 randomized splits, e.g., a 80/20 train/test split repeated 5 times, with non-overlapping splits
# Here we instantiate the KFold CV-splitter object so that we can define the random seed for replicability
y_pred = cross_val_predict(pipe, x, y, cv=KFold(n_splits=5, shuffle=True, random_state=42))
title = "Cross-validation, 5 randomized splits \n" r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred), r2_score(y, y_pred))
plotyy(y, y_pred, color=eis_data.seriesIdx, title=title)<jupyter_output><empty_output><jupyter_text>Note that while we are plotting the entire data set, every one of these predictions is a 'test' prediction. Again, with splits randomized by data point and no domain context (knowledge of which cell each data point is from), it looks like we have a very good model. We can test how well the model performns when being split on a cell basis, i.e., our seriesIdx variable, using functionality already built in.<jupyter_code>from sklearn.model_selection import LeaveOneGroupOut

# Instatiate the pipeline object
pipe = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('regressor', LinearRegression())
    ]
)
# Also make a dummy pipeline
dummy_pipe = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('regressor', DummyRegressor())
    ]
)
# Run cross-validation by group, using a leave-on
y_pred = cross_val_predict(pipe, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())
title = "Cross-validation, leave-one-out by group \n" r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred), r2_score(y, y_pred))
plotyy(y, y_pred, color=eis_data.seriesIdx, title=title)
# Dummy results
y_pred_dummy = cross_val_predict(dummy_pipe, x, y, groups = eis_data.seriesIdx, cv=LeaveOneGroupOut())
title = "Dummy CV, leave-one-out by group \n" r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_dummy), r2_score(y, y_pred_dummy))
plotyy(y, y_pred_dummy, color=eis_data.seriesIdx, title=title)<jupyter_output><empty_output><jupyter_text>This is dramatically worse. Again, the dummy model, which simply uses the mean of the training set to make predictions, is actually more accurate on average. In summary, so far:- Linear model error when training on the entire data set is extremely low (MAE = 0.36 mAh)- Linear model error during leave-one-out cross-validation by group is ...    - ...several orders of magnitude worse than training (MAE = 4.27 mAh)    - ...worse than using a dummy model (dummy MAE = 3.85 mAh)These three observations are classic signs that a model is overfit to the data. Overfitting results in extremely poor performance on test data. There are many ways to improve our model performance on test data (also called 'regularizing' the model). Two methods we will explore here are:1. choosing a type of estimator (regressor object) that is less prone to overfitting than a linear model, and2. feature engineering Model regularization: using a different regressorSimply changing the type of regression model used can dramatically impact results. Here, we will compare three options:- Linear regression- Gaussian process regression- XGBoost regressionThese three models have substantially different approaches, and should result in different final outcomes. Dramatically simplifying these models: linear models use least-squares minimization to optimize a bias and slopes for each feature; GPR models utilize a kernel function to generate a prior and use log-likelihood minimization to optimize the posterior distribution; XGBoost is a tree based model, using randomization to try and prevent overfitting.<jupyter_code>from sklearn.gaussian_process import GaussianProcessRegressor
from xgboost import XGBRegressor

# Always compare to a dummy
dummy_pipe = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('regressor', DummyRegressor())
    ]
)
# Linear regression
pipe_linear = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('regressor', LinearRegression())
    ]
)
# Gaussian process regression
pipe_gpr = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('regressor', GaussianProcessRegressor())
    ]
)
# XGBoost regressor
pipe_xgb = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('regressor', XGBRegressor())
    ]
)

# Dummy
y_pred_dummy = cross_val_predict(dummy_pipe, x, y, groups = eis_data.seriesIdx, cv=LeaveOneGroupOut())
mae_dummy = mean_absolute_error(y, y_pred_dummy)
# Linear
y_pred_linear = cross_val_predict(pipe_linear, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())
mae_linear = mean_absolute_error(y, y_pred_linear)
# Gaussian process
y_pred_gpr = cross_val_predict(pipe_gpr, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())
mae_gpr = mean_absolute_error(y, y_pred_gpr)
# XGBoost
y_pred_xgb = cross_val_predict(pipe_xgb, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())
mae_xgb = mean_absolute_error(y, y_pred_xgb)
regressors = ['Dummy','Linear','GPR','XGBoost']
maes = [mae_dummy, mae_linear, mae_gpr, mae_xgb]
plt.bar(regressors, maes)
plt.ylabel(r'MAE$_{Test}$ (mAh)')
plt.show()<jupyter_output><empty_output><jupyter_text>Only XGBoost is better than the dummy regressor. Visualize the results.<jupyter_code>title = "XGBoost, CV leave-one-out by group \n" r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_xgb), r2_score(y, y_pred_xgb))
plotyy(y, y_pred_xgb, color=eis_data.seriesIdx, title=title)<jupyter_output><empty_output><jupyter_text>This is quite a lot better than the linear model, the error is ~10% of the average BOL capacity (35 mAh), and all the trends look correct.  Model regularization: feature engineeringFeature engineering is a broad term with covering many possible approaches for handling features. Generalizing a bit, there are perhaps three overarching types of feature engineering:- Feature generation (few features -> many)    - Ex.: X__i,j = X_i * X_j for all combinations of features i, j- Feature selection (many features -> few)    - Ex.: Use impedance data from only 100 Hz- Feature extraction (many interrelated features -> few independent features)    - Ex.: Fit ECM to impedance data, use ECM parameters as featuresHere, we will study two of the options: feature selection and feature extraction. First, feature selection. We can simply consider a feature selection algorithm that uses only the values from a single frequency. We can automate the extraction of this frequency by defining a class, and passing it to the Pipeline. This might seem like a fair bit of work for a simple operation, but the benefit here is that this is now fully automated, and can happen at any point in an arbitrarily complex sequence of feature transformation operations.<jupyter_code>from sklearn.base import BaseEstimator, TransformerMixin
# Inherent sklearn classes for transformers, this gives us a bunch of methods for free, see sklearn documentation
class SingleFrequencySelector(BaseEstimator, TransformerMixin): 
    def __init__(self, freq, index=0):
        self.freq = freq
        self.index = index
    
    def fit(self, X, y=None):
        return self
    
    def transform(self, X, y=None):
        X_ = X.copy()
        n = self.freq.size
        index = [self.index, self.index+n, self.index+n*2, self.index+n*3]
        return X_[:,index]<jupyter_output><empty_output><jupyter_text>Test just picking the first frequency.<jupyter_code># Linear regression, single frequency
index_freq = 0
pipe_lin_1freq = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('selector', SingleFrequencySelector(freq=freq, index=index_freq)),
        ('regressor', LinearRegression())
    ]
)
# Linear
y_pred_lin_1freq = cross_val_predict(pipe_lin_1freq, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())
title = "Linear (Frequency: %0.2g Hz) \n" % (freq[index_freq]) + "CV leave-one-out by group \n" + r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_lin_1freq), r2_score(y, y_pred_lin_1freq))
plotyy(y, y_pred_lin_1freq, color=eis_data.seriesIdx, title=title)<jupyter_output><empty_output><jupyter_text>This wasn't a good frequency to pick - we need to find out if picking a better frequency results in a better model. But we don't want to manually try all 60 frequencies.One of the benefits of automation is that we can take advantage of sklearn functionality to automatically grid-search across all possible frequencies.<jupyter_code>from sklearn.model_selection import GridSearchCV
# Linear regression, single frequency
pipe_lin_1freq = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('selector', SingleFrequencySelector(freq=freq)),
        ('regressor', LinearRegression())
    ]
)
# Parameters of any pipeline step can be swept using the name of the step, a dunder '__', and the name of the parameter
param_grid = {
    "selector__index": np.arange(0, freq.size)
}
search = GridSearchCV(pipe_lin_1freq,
    param_grid=param_grid,
    cv=LeaveOneGroupOut(),
    scoring='neg_mean_absolute_error'
)
search.fit(x, y, groups=eis_data.seriesIdx)

# We can make model predictions using the best result exactly the same as before
y_pred_lin_1freq = cross_val_predict(search.best_estimator_, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())
title = "Linear (Frequency: %0.2g Hz) \n" % (freq[search.best_params_['selector__index']]) + "CV leave-one-out by group \n" + r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_lin_1freq), r2_score(y, y_pred_lin_1freq))
plotyy(y, y_pred_lin_1freq, color=eis_data.seriesIdx, title=title)<jupyter_output><empty_output><jupyter_text>This is better than the XGBoost model, which has a much more complicated model architecture. All we had to do was find a good frequency. But we're still making some extreme predictions near end-of-life.We can investigate the MAE versus frequency to help see how the frequency impacts model performance.<jupyter_code>plt.semilogx(freq, -1*search.cv_results_['mean_test_score'], '-k')
plt.axvline(freq[search.best_params_['selector__index']], color='r', linestyle='--')
plt.xlabel('Frequency (Hz)')
plt.ylabel(r'MAE$_{Cross-validation}$')
plt.legend(['All models', 'Best'])
plt.title('Linear models, single frequency')<jupyter_output><empty_output><jupyter_text>Note that we are basically constructing our own feature selection method here, and optimizing the hyperparameters of the feature selection algorithm with an exhaustive grid search. For longer searches, sklearn has other tools such as `RandomizedSearchCV`, `HalvingGridSearchCV`, and `HalvingRandomSearchCV`.To get an idea of how a model is making a prediction, there are many tools we can use to iterrogate a model. Some popular examples include SHAP for estimating contributions of each feature, or partial dependence plots. Note that both of these methods are model agnostic - there are also model specific methods. Here we show partial dependence. Partial dependence shows how the average output value varies when the value of a specific feature is set to a single value across the entire data set and varied across the observed range. We also plot the predicted value from every single row of the data set to show how the results vary independently of the feature being probed by the partial dependence calculation.<jupyter_code>from sklearn.inspection import PartialDependenceDisplay
index_best = search.best_params_['selector__index']
n = freq.size
index_best = [index_best, index_best+n, index_best+n*2, index_best+n*3]
fig, ax = plt.subplots(1, 4, figsize=(12,3))
PartialDependenceDisplay.from_estimator(search.best_estimator_,
    x, features=index_best, kind='both',
    ice_lines_kw = {'linestyle': '-', 'color': 'k', 'linewidth': 1},
    pd_line_kw = {'linestyle': '--', 'color': 'r', 'linewidth': 2},
    ax=ax,
)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>From above, the Z_Real and |Z| features seem by far the most important, but note the y-scale: the impact is huge, causing the expected value of the output to be extreme. Also, our two most important features are cancelling each other out - this implies our model is extremely sensitive, because if we get some values of those two features that are not well covered in our training data, the model prediction may be extremely bad. This is likely what is causing the very poor predictions at low capacity.Partial dependence can also be done on two features simultaneously, performing a 2D grid search. We can try this on our two most important features to see what happens.<jupyter_code>plt.set_cmap('cividis')
PartialDependenceDisplay.from_estimator(search.best_estimator_,
    x, features=[(index_best[0], index_best[2])],
)
# Scatter of feature values from the actual data
feature = x.columns[[index_best[0], index_best[2]]]
n_series = np.unique(eis_data.seriesIdx).size
cmap = plt.get_cmap('tab10', n_series)
sc = plt.scatter(x[feature[0]].values, x[feature[1]].values, c=eis_data.seriesIdx, cmap=cmap)
cbar = plt.colorbar(sc)
tick_locs = (np.arange(n_series) + 1.5)*(n_series-1)/n_series
cbar.set_ticks(tick_locs)
cbar.set_ticklabels(np.arange(n_series)+1)<jupyter_output><empty_output><jupyter_text>Points corresponding to the training data have been added to the partial dependence plot. We can see that these two features are highly correlated, and as noted, they have opposing slopes (resulting in the diagonal trend for the partial dependence). Basically, this model is a see-saw, with extremely heavy weights on either side, which is currently precariously balanced to give accurate predictions in an extremely narrow regime - if the value of either parameter varies even a tiny bit from the distribution of the training data, this model will make an extremely poor prediction.We can compare this partial dependence plot to one calculated using an XGBoost estimator to see how different the model structure is. Let's just use the same frequencies, though XGBoost may have a different set of optimal frequencies compared to a linear model.<jupyter_code># XGBoost regression, single frequency
index_freq = search.best_params_['selector__index']
pipe_xgb_1freq = Pipeline(
    [
        ('scaler', StandardScaler()),
        ('selector', SingleFrequencySelector(freq=freq, index=index_freq)),
        ('regressor', XGBRegressor())
    ]
)
# Linear
y_pred_xgb_1freq = cross_val_predict(pipe_xgb_1freq, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())
title = "XGBoost (Frequency: %0.2g Hz) \n" % (freq[index_freq]) + "CV leave-one-out by group \n" + r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_xgb_1freq), r2_score(y, y_pred_xgb_1freq))
plotyy(y, y_pred_xgb_1freq, color=eis_data.seriesIdx, title=title)
from sklearn.inspection import PartialDependenceDisplay
index_best = search.best_params_['selector__index']
n = freq.size
# Need to fit the model on all data first
pipe_xgb_1freq.fit(x, y)
# Then plot partial dependencies
index_best = [index_best, index_best+n, index_best+n*2, index_best+n*3]
fig, ax = plt.subplots(1, 4, figsize=(12,3))
PartialDependenceDisplay.from_estimator(pipe_xgb_1freq,
    x, features=index_best, kind='both',
    ice_lines_kw = {'linestyle': '-', 'color': 'k', 'linewidth': 1},
    pd_line_kw = {'linestyle': '--', 'color': 'r', 'linewidth': 2},
    ax=ax,
)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>The relationships are obviously non-linear now, and the expected values (y-axis) have become way more reasonable. The second and fourth features seem like the most important.<jupyter_code>plt.set_cmap('cividis')
PartialDependenceDisplay.from_estimator(pipe_xgb_1freq,
    x, features=[(index_best[1], index_best[3])],
)
# Scatter of feature values from the actual data
feature = x.columns[[index_best[1], index_best[3]]]
n_series = np.unique(eis_data.seriesIdx).size
cmap = plt.get_cmap('tab10', n_series)
sc = plt.scatter(x[feature[0]].values, x[feature[1]].values, c=eis_data.seriesIdx, cmap=cmap)
cbar = plt.colorbar(sc)
tick_locs = (np.arange(n_series) + 1.5)*(n_series-1)/n_series
cbar.set_ticks(tick_locs)
cbar.set_ticklabels(np.arange(n_series)+1)<jupyter_output><empty_output><jupyter_text>Again, the relationships are highly non-linear. It seems like there's two separated 'low capacity' regions, which may correspond to cells showing different behaviors in the training set - there seem to be two 'groups' of cells, with data from most cells near the top of the plot, and data from 2 cells near the bottom. But, unlike the linear model, the predictions are much more consistent, not varying wildly outside of the expected range, indicating this model may be more robust to outliers or use in the real-world than the linear model. It's unlikely that this is the absolute best we can do. We can try other models (support vector regression, XGBoost, GPR), extracting multiple frequencies instead of 1, or trying to use domain knowledge (such as pulling interesting points of the Nyquist plot). For example, sklearn has built-in feature selection methods, such as `SelectKBest`, which only keeps the *k* features that are most highly scored versus the target variable, using some scoring function. Compare current results:<jupyter_code>from sklearn.metrics import max_error

# x labels
regressors = ['Dummy',
    'Linear',"Linear \n (%0.2g Hz)" % (freq[search.best_params_['selector__index']]),
    'XGBoost', "XGBoost \n (%0.2g Hz)" % (freq[search.best_params_['selector__index']])]
# mean absolute errors
mae_lin_1freq = mean_absolute_error(y, y_pred_lin_1freq)
mae_xgb_1freq = mean_absolute_error(y, y_pred_xgb_1freq)
maes = [mae_dummy, mae_linear, mae_lin_1freq, mae_xgb, mae_xgb_1freq]
# max absolute errors
max_errors = [
    max_error(y, y_pred_dummy),
    max_error(y, y_pred_linear),
    max_error(y, y_pred_lin_1freq),
    max_error(y, y_pred_xgb),
    max_error(y, y_pred_xgb_1freq),
]

# Plot
fig, ax1 = plt.subplots()

b1 = ax1.bar(regressors, maes, color='tab:blue', edgecolor='k', label='MAE', width=-0.3, align='edge')
ax1.bar_label(b1, fmt='%.1f')
ax1.set_ylabel(r'MAE$_{Test}$ (mAh)', color='tab:blue')

ax2 = ax1.twinx()
b2 = plt.bar(regressors, max_errors, color='tab:orange', edgecolor='k', label='MaxAE', width=0.3, align='edge')
ax2.bar_label(b2, fmt='%.1f')
ax2.set_ylabel(r'MaxAE$_{Test}$ (mAh)', color='tab:orange')

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Even though the average error of the linear single frequency model is the best, the maximum error is very large, indicating the high sensitivity of the model that we found using partial dependence. The best result is the model with both low average error and low maximum error, in this case, the XGBoost model. Given that choosing a best model often balances multiple metrics and domain considerations (computation time for training or making predictions, memory requirements), it is not always obvious what the 'best model' is.  Summary- Process data into a format for easy input into machine-learning libraries and/or visualization tools- Visualize your target variable(s) and feature(s) as best you can    - Use domain knowledge to help interpret the complexity of your data and the modeling challenge    - Consider physical relationships between variables that could guide you modeling effort    - Be careful of data leakage (spurious correlations, i.e., is your model learning physical relationships or learning how you set up your experiment/processed your data?)    - Outliers can substantially impact modeling results, and are often not obvious- Split you data into training, test, and validation splits that address a real-world relevant challenge    - Random train/test splits are often not challenging or interesting    - If there's not enough data to split into three sets, consider using cross-validation or bootstrap resampling on your training set- Try many possible approaches before assuming any one clever idea is best    - Use automation to optimize hyperparameters of any particular approach- Interrogate models by visualizing results in domain context and using model interrogation methods to try and develop an intuition for model behavior- Use multiple error metrics when deciding on a 'best' model  ChallengeFind a model with lower error than the linear, single frequency model above. Some ideas:- Search for the best combination of two frequencies- Use a more complex feature selection strategy     - Ex., using feature importance from a `RandomForestRegressor`, see `SelectFromModel`- Use a different estimator that better handles non-linear relationships, combined with feature selection approaches    - Ex., `XBGRegressor`, `SVR`, `GaussianProcessRegressor`- Extract features using domain knowledge    - Ex., track minima or maxima from Nyquist plot<jupyter_code># Search for the best combination of two frequencies<jupyter_output><empty_output><jupyter_text>Click for code hints for above```class DoubleFrequencySelector(BaseEstimator, TransformerMixin):     def __init__(self, freq, index_0=0, index_1=0):        self.freq = freq        self.index_0 = index_0        self.index_1 = index_1        def fit(self, X, y=None):        return self        def transform(self, X, y=None):        X_ = X.copy()        n = self.freq.size        index_0 = [self.index_0, self.index_0+n, self.index_0+n*2, self.index_0+n*3]        index_1 = [self.index_1, self.index_1+n, self.index_1+n*2, self.index_1+n*3]        index = np.sort(np.concatenate((index_0,index_1)))        return X_[:,index]pipe_lin_2freq = Pipeline(    [        ('scaler', StandardScaler()),        ('selector', DoubleFrequencySelector(freq=freq)),        ('regressor', LinearRegression())    ])param_grid = {    "selector__index_0": np.arange(0, freq.size, step=3),    "selector__index_1": np.arange(1, freq.size, step=3),  avoid redundant features}search = GridSearchCV(pipe_lin_2freq,    param_grid=param_grid,    cv=LeaveOneGroupOut(),    scoring='neg_mean_absolute_error')search.fit(x, y, groups=eis_data.seriesIdx)y_pred_lin_2freq = cross_val_predict(search.best_estimator_, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())title = "Linear (%0.2g Hz, %0.2g Hz) \n" % (freq[search.best_params_['selector__index_0']], freq[search.best_params_['selector__index_1']]) + "CV leave-one-out by group \n" + r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_lin_2freq), r2_score(y, y_pred_lin_2freq))plotyy(y, y_pred_lin_2freq, color=eis_data.seriesIdx, title=title)```<jupyter_code># Use a more complex feature selection strategy, ex., SelectFromModel<jupyter_output><empty_output><jupyter_text>Click for code hints for above```from sklearn.ensemble import ExtraTreesRegressorfrom sklearn.feature_selection import SelectFromModel Linear regression, select from modelpipe_lin_treefeatures = Pipeline(    [        ('scaler', StandardScaler()),        ('selector', SelectFromModel(ExtraTreesRegressor())),        ('regressor', LinearRegression())    ])y_pred_lin_treefeatures = cross_val_predict(pipe_lin_treefeatures, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())title = "Linear (tree features) \n CV leave-one-out by group \n" + r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_lin_treefeatures), r2_score(y, y_pred_lin_treefeatures))plotyy(y, y_pred_lin_treefeatures, color=eis_data.seriesIdx, title=title)```<jupyter_code># Try a different estimator with feature selection, ex., SVR with single feature selection<jupyter_output><empty_output><jupyter_text>Click for code hints for above```from sklearn.svm import SVRpipe_svr_1freq = Pipeline(    [        ('scaler', StandardScaler()),        ('selector', SingleFrequencySelector(freq=freq)),        ('regressor', SVR())    ])param_grid = {    "selector__index": np.arange(0, freq.size, 3),    "regressor__C": np.logspace(-4, 1, 10),}search = GridSearchCV(    pipe_svr_1freq,    param_grid=param_grid,    cv=LeaveOneGroupOut(),    scoring='neg_mean_absolute_error')search.fit(x, y, groups=eis_data.seriesIdx)y_pred_svr_1freq = cross_val_predict(search.best_estimator_, x, y, groups=eis_data.seriesIdx, cv=LeaveOneGroupOut())title = "SVR (Frequency: %0.2g Hz) \n" % (freq[search.best_params_['selector__index']]) + "CV leave-one-out by group \n" + r"MAE: %0.2f, R$^2$: %0.2f" % (mean_absolute_error(y, y_pred_svr_1freq), r2_score(y, y_pred_svr_1freq))plotyy(y, y_pred_svr_1freq, color=eis_data.seriesIdx, title=title)```<jupyter_code># Extract features using domain knowledge, ex., track the valley in Z_imaginary at low frequency<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09413395165374057, 'fraction_numerical': 0.010527437968071214, 'mean_word_length': 5.541492329149233, 'pattern_counts': {'":': 5, '<': 123, '<?xml version=': 0, '>': 126, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '615051', 'license': 'MIT', 'n_tokens_mistral': 11777, 'n_tokens_neox': 11298, 'n_words': 4057, 'path': 'predict_capacity_from_eis/ml_predicting_capacity_with_eis.ipynb', 'repo_name': 'battery-data-commons/mrs-sp22-tutorial', 'size': '37540'}