<jupyter_start><jupyter_text>Train a Model to Predict Formation Energy using the OQMDThis notebook recreates a 2016 paper by [Ward et al.](https://www.nature.com/articles/npjcompumats201628) on predicting the formation enthalpy of materials based on their composition. We will use the [Materials Data Facility](http://materialsdatafacility.org) to retrieve a training set from the the [OQMD](http://oqmd.org), compute features based on the composition of each entry, and then train a random forest model.This example was last updated on 11/15/18 for Matminer v.0.4.5<jupyter_code>%matplotlib inline
from matminer.data_retrieval import retrieve_MDF
from matminer.featurizers.base import MultipleFeaturizer
from matminer.featurizers import composition as cf
from matminer.featurizers.conversions import StrToComposition
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm
import numpy as np
import pandas as pd
import pickle as pkl
from sklearn import metrics
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score, cross_val_predict, GridSearchCV, ShuffleSplit, KFold<jupyter_output><empty_output><jupyter_text>Settings to change<jupyter_code>quick_demo = True # Whether to run an faster version of this demo. 
# The full OQMD model takes about a hour to test and ~8GB of RAM<jupyter_output><empty_output><jupyter_text>Load Training SetWard _et al._ trained their machine learning models on the formation enthalpies of crystalline compounds form the [OQMD](oqmd.org). Here, we extract the data using the copy of the OQMD available through the MDF  Download the DataWe first create a `Forge` instance, which simplifies performing search queries against the MDF. The first step is to create a tool for reading from the MDF's search index.<jupyter_code>mdf = retrieve_MDF.MDFDataRetrieval(anonymous=True)<jupyter_output><empty_output><jupyter_text>Then, we assemble a query that gets only the converged static calculations from the OQMD.<jupyter_code>query_string = 'mdf.source_name:oqmd AND (oqmd.configuration:static OR '\
    'oqmd.configuration:standard) AND dft.converged:True'
if quick_demo:
    query_string += " AND mdf.scroll_id:<10000"
data = mdf.get_data(query_string, unwind_arrays=False)<jupyter_output><empty_output><jupyter_text>This tool creates a DataFrame object with the metadata for each entry in the OQMD<jupyter_code>data.head(2)<jupyter_output><empty_output><jupyter_text>We only need two columns: `delta_e` and `material.composition`<jupyter_code>data = data[['oqmd.delta_e.value', 'material.composition']]<jupyter_output><empty_output><jupyter_text>Renaming the columns to make the rest of the code more succinct<jupyter_code>data = data.rename(columns={'oqmd.delta_e.value': 'delta_e', 'material.composition':'composition'})<jupyter_output><empty_output><jupyter_text>Compile the Training SetOur next step is to get only the lowest-energy entry for each composition.<jupyter_code>data = StrToComposition(target_col_id='composition_obj').featurize_dataframe(data, 'composition')<jupyter_output><empty_output><jupyter_text>Create shortcuts for our input and output columns Remove compounds w/o a `delta_e` measurement.<jupyter_code>for k in ['delta_e']:
    data[k] = pd.to_numeric(data[k])
original_count = len(data)
data = data[~ data['delta_e'].isnull()]
print('Removed %d/%d entries'%(original_count - len(data), original_count))<jupyter_output>Removed 409/4849 entries<jupyter_text>Get only the groundstate and each composition<jupyter_code>%%time
original_count = len(data)
data['composition'] = data['composition_obj'].apply(lambda x: x.reduced_formula)
data.sort_values('delta_e', ascending=True, inplace=True)
data.drop_duplicates('composition', keep='first', inplace=True)
print('Removed %d/%d entries'%(original_count - len(data), original_count))<jupyter_output>Removed 24/4440 entries
CPU times: user 643 ms, sys: 8.83 ms, total: 652 ms
Wall time: 742 ms<jupyter_text>Remove outliers<jupyter_code>original_count = len(data)
data = data[np.logical_and(data['delta_e'] >= -20, data['delta_e'] <= 5)]
print('Removed %d/%d entries'%(original_count - len(data), original_count))<jupyter_output>Removed 1/4416 entries<jupyter_text>Build an ML modelIn this part of the notebook, we build a ML model using [scikit-learn](http://scikit-learn.org/stable/) and evaluate its performance using cross-validation.  Part 1: Compute RepresentationThe first step in building a ML model is to convert the raw materials data (here: the composition) into the required input for an ML model: a finite list of quantitative attributes. In this example, we use the "general-purpose" attributes of [Ward *et al* 2016](https://www.nature.com/articles/npjcompumats201628).<jupyter_code>feature_calculators = MultipleFeaturizer([cf.Stoichiometry(), cf.ElementProperty.from_preset("magpie"),
                                          cf.ValenceOrbital(props=['avg']), cf.IonProperty(fast=True)])<jupyter_output><empty_output><jupyter_text>Get the feature names<jupyter_code>feature_labels = feature_calculators.feature_labels()<jupyter_output><empty_output><jupyter_text>Compute the features<jupyter_code>%%time
data = feature_calculators.featurize_dataframe(data, col_id='composition_obj');
print('Generated %d features'%len(feature_labels))
print('Training set size:', 'x'.join([str(x) for x in data[feature_labels].shape]))<jupyter_output>Generated 145 features
Training set size: 4415x145<jupyter_text>Remove entries with `NaN` or `infinite` features<jupyter_code>original_count = len(data)
data = data[~ data[feature_labels].isnull().any(axis=1)]
print('Removed %d/%d entries'%(original_count - len(data), original_count))<jupyter_output>Removed 0/4415 entries<jupyter_text>Part 2: Tuning Model Hyper-parametersFor brevity, we will only consider one ML algorithm in this example: [random forest](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.htmlsklearn.ensemble.RandomForestRegressor). The "random forest" algorithm works by training many different decision tree models, where each is trained on a different subset of the dataset . Here, we tune one of the major parameters of the algorithm: the number features considered at each split in each decision tree<jupyter_code>model = GridSearchCV(RandomForestRegressor(n_estimators=20 if quick_demo else 150, n_jobs=-1),
                     param_grid=dict(max_features=range(8,15)),
                     scoring='neg_mean_squared_error',cv=ShuffleSplit(n_splits=1, test_size=0.1))
model.fit(data[feature_labels], data['delta_e'])<jupyter_output><empty_output><jupyter_text>Plot the tuning results. This shows the CV score as a function of the parameter we tuned "max features"<jupyter_code>model.best_score_
fig, ax = plt.subplots()

# Plot the score as a function of alpha
ax.scatter(model.cv_results_['param_max_features'].data,
          np.sqrt(-1 * model.cv_results_['mean_test_score']))
ax.scatter([model.best_params_['max_features']], np.sqrt([-1*model.best_score_]), marker='o', color='r', s=40)
ax.set_xlabel('Max. Features')
ax.set_ylabel('RMSE (eV/atom)')<jupyter_output><empty_output><jupyter_text>Save our best model<jupyter_code>model = model.best_estimator_<jupyter_output><empty_output><jupyter_text>Part 3: Cross-validation TestQuantify the performance of this model using 10-fold cross-validation<jupyter_code>cv_prediction = cross_val_predict(model, data[feature_labels], data['delta_e'], cv=KFold(10, shuffle=True))<jupyter_output><empty_output><jupyter_text>Compute aggregate statistics<jupyter_code>for scorer in ['r2_score', 'mean_absolute_error', 'mean_squared_error']:
    score = getattr(metrics,scorer)(data['delta_e'], cv_prediction)
    print(scorer, score)
model<jupyter_output><empty_output><jupyter_text>Plot the individual predictions<jupyter_code>fig, ax = plt.subplots()

ax.hist2d(pd.to_numeric(data['delta_e']), cv_prediction, norm=LogNorm(), bins=64, cmap='Blues', alpha=0.9)

ax.set_xlim(ax.get_ylim())
ax.set_ylim(ax.get_xlim())

mae = metrics.mean_absolute_error(data['delta_e'], cv_prediction)
r2 = metrics.r2_score(data['delta_e'], cv_prediction)
ax.text(0.5, 0.1, 'MAE: {:.2f} eV/atom\n$R^2$:  {:.2f}'.format(mae, r2),
        transform=ax.transAxes,
       bbox={'facecolor': 'w', 'edgecolor': 'k'})

ax.plot(ax.get_xlim(), ax.get_xlim(), 'k--')

ax.set_xlabel('DFT $\Delta H_f$ (eV/atom)')
ax.set_ylabel('ML $\Delta H_f$ (eV/atom)')

fig.set_size_inches(3, 3)
fig.tight_layout()
fig.savefig('oqmd_cv.png', dpi=320)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11214953271028037, 'fraction_numerical': 0.015070093457943926, 'mean_word_length': 7.0688030160226205, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 81, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '284626', 'license': 'BSD-3-Clause-LBNL', 'n_tokens_mistral': 2949, 'n_tokens_neox': 2800, 'n_words': 697, 'path': 'matminer_examples/machine_learning-nb/formation_e.ipynb', 'repo_name': 'flaviu-gostin/matminer_examples', 'size': '8560'}