<jupyter_start><jupyter_text>Digging Into Deep Time and Deep Cover<NAME> , <NAME>, <NAME> and <NAME>; CSIRO Mineral Resources Contents| [**Abstract**](./00_overview.ipynb) | **Introduction**                                                    | [**Examples**](./00_overview.ipynbExamples)            | [**Tools**](./00_overview.ipynbLeveraging-the-Scientific-Python-Ecosystem) | |:-----|:-----|:-----|:-----||  | [Minerals Exploration](./00_overview.ipynbAn-Evolving-Role-of-Geochemistry-in-Mineral-Exploration)  | [Classification](./011_classification.ipynb) |  ||  | [Data Driven Geochem](./00_overview.ipynbData-Driven-Geochemistry) | [Data Visualization](./012_datavis.ipynb) |  |  Classification: Towards Multivariate ClassifiersIn geochemistry, classification is for historical reasons often treated as a graphical or two-dimensional problem for which the output includes a graphically-represented classification or discrimination diagram (driven by a need to visualize the data; e.g. [Pearce, 1973]). And while many of these diagrams have been practically useful for 'binning' and describing rocks, the divisions we use don't necessarily make sense in the natural world (e.g. consider the Total Alkali - Silica diagram). Modern statistical and machine learning techniques can readily be applied to make better use of the data dimensionality and move beyond the human-centric limitations. In this case a data driven approach instead considers divisions according to the natural clustering of data, where they exist.**Key Messages**:* Use the dimensionality you have* Particularly if you're using major element data, use compositional transforms* Where available, use metrics to understand estimated uncertainty and confidence* Try to visualize and understand what's happening in lower dimensional spaces Searching for Enriched BasaltsThis example is motivated by one of the research questions for my postdoc: where and when do we find the earliest examples of enriched basalts (analogous to modern day OIB)? Large magmatic nickel deposits are assocaited with enriched basaltic suites, and if we can generate spatio-temporal constraints on their genesis (particularly in geological time), we might be able to reduce the nickel exploration search space considerably. One approach to this question is to formulate it as a classification problem. Whole rock geochemistry is widely used to constrain the tectonic context in which magmatic rocks have been formed, with varying levels of success. Tectonic discrimination using whole-rock geochemical data is complicated by long geological histories, alteration, biased preservation and secular change. Notably, there is a lack of confidently identified oceanic crust in the Archean, suggesting finding OIB-like rocks is likely even more challenging. Most traditional tectonic discrimination diagrams use geochemical proxies to encode information about geological processes or (e.g. the "Pearce diagrams" which use Th/Yb, Nb/Yb and TiO2/Yb as proxies for subduction influence, source enrichment and depth of mantle melting respectively; [Pearce, 2008]), while others are designed to provide optimal discrimination for a specific scenario. While a wide variety of geochemical proxies can be used with relative confidence to discriminate between modern rocks, investigating primary geochemical contrast in ancient rocks (especially into the Archean) is typically considerd to require restriction of variables considered to those which are least effected by geological modification, i.e. immobile elements such as the High Field Strength Elements (HFSE) and Rare Earth Elements (REE). Overall the limited perspective of graphically-focused methods can lead to variable and uncertain classifications depending on the elements chosen, especially when applied into deep time. The example below goes through some of the first steps of exploring how we might be able to classify Paleoproteroizic and Archean basalts based on their geochemistry with a machine learning approach which capitalizes on the dimensionality of our datasets. But before we start digging into deep time (also see the [note on dealing with space and time](./00_overview.ipynbOther-Opportunities-in-Data-Integration-and-Exploration) in the overview notebook), we first consider how well we can distinguish modern ocean island basalts. More so than the results of the prediction, we'd like to demonstrate and explore part of the process. Machine Tectonic DiscriminationBelow we'll demonstrate a machine learning approach to classifying basalt geochemistry based on tectonic setting. We'll adapt a version of a support vector classifier, similar to those previously used for this problem ([Petrelli, 2016]; [Ueki, 2018]), and re-use the same 29-dimension tectonic discrimination dataset (just over 2000 samples taken from global repositories, each with majors, minors, traces and isotopes) from Ueki et al. (2018), modified to record the tectonic setting of each basalt in a 'Class' column.First let's examine a well-used trace element discrimination diagram (the Th-Nb-Yb diagram; [Pearce, 2008]) to get some insight what some of the key issues with classification using only a few elements or dimension. This is one of number of classification and tectonic diagrams have been developed which use specific geochemical proxies to link rock and mineral geochemistry to geological processes and reservoirs (here Nb/Yb and Th/Yb which distinguish relative enrichment of the mantle and arc-related processes). As the diagram presents only two dimensions (or three, depending who you talk to), the data presented show only part of the whole picture. This reduced dimensionality means that when you compare rocks from a variety of settings, you're likely to see large degrees of overlap, as found by [Li et al. (2015)], and shown below. While the diagram is not particularly useful as a classification tool, the proxies used (ratios or otherwise) still provide useful 'features' which can be added to datasets before constructing a multivariate classifier.| | ||--|--|| Figure from [Smithies (2018)] illustrating the main features of the Th-Yb-Nb diagram [Pearce (2008)], including the mantle or 'MORB-OIB array' and the 'modern arc array'. | Figure from [Li et al. (2015)], illustrating the significant degree of overlap in these two dimensions between different tectonic settings renders this approach futile for generalized discrimination.  |Below we've loaded the dataset from [Ueki, 2018] and added the Nb/Yb and Th/Yb, and contoured data density for basalts from three different tectonic settings (at 30th, 50th and 80th percentiles). Notably, these data density contours largely line up with what we expect from the Pearce diagram.[Pearce, 1973]: https://doi.org/10.1016/0012-821X(73)90129-5 "<NAME>., <NAME>., 1973. Tectonic setting of basic volcanic rocks determined using trace element analyses. Earth and Planetary Science Letters 19, 290–300. doi: 10.1016/0012-821X(73)90129-5"[Li et al. (2015)]: https://doi.org/10.1016/j.lithos.2015.06.022 "<NAME>., <NAME>., <NAME>., <NAME>., 2015. Trace element indiscrimination diagrams. Lithos 232, 76–83. "[Smithies (2018)]: https://doi.org/10.1016/j.epsl.2018.01.034 "<NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., <NAME>., 2018. Two distinct origins for Archean greenstone belts. Earth and Planetary Science Letters 487, 106–116."[Pearce (2008)]: https://doi.org/10.1016/j.lithos.2007.06.016 "<NAME>., 2008. Geochemical fingerprinting of oceanic basalts with applications to ophiolite classification and the search for Archean oceanic crust. Lithos 100, 14–48."[Pearce, 2008]: https://doi.org/10.1016/j.lithos.2007.06.016 "<NAME>., 2008. Geochemical fingerprinting of oceanic basalts with applications to ophiolite classification and the search for Archean oceanic crust. Lithos 100, 14–48."[Petrelli, 2016]: https://doi.org/10.1007/s00410-016-1292-2 "<NAME>., <NAME>., 2016. Solving petrological problems through machine learning: the study case of tectonic discrimination using geochemical and isotopic data. Contrib Mineral Petrol 171, 81."[Ueki, 2018]: https://doi.org/10.1029/2017GC007401 "<NAME>., <NAME>., <NAME>., 2018. Geochemical Discrimination and Characteristics of Magmatic Tectonic Settings: A Machine-Learning-Based Approach. Geochemistry, Geophysics, Geosystems 19, 1327–1347."<jupyter_code>import pandas as pd
import pyrolite.geochem

df = pd.read_csv("https://storage.googleapis.com/aegc2019/ueki2018.csv")
df.pyrochem.add_ratio("Th/Yb")
df.pyrochem.add_ratio("Nb/Yb");
import pyrolite.plot
import matplotlib.pyplot as plt
from pyrolite.plot.templates.pearce import pearceThNbYb
from pyrolite.util.plot.legend import proxy_line

fig, ax = plt.subplots(1, figsize=(8, 6))
pearceThNbYb(ax)

proxies = []
labels = []
for setting, ID, cmap in [("MORB", "MOR", "Blues"),
                          ("Continental Arcs", "CA","Reds"),
                          ("Ocean Islands", "OI", "Greens")]:
    df.loc[(df.Class == ID), ["Nb/Yb", "Th/Yb"]].pyroplot.density(
        ax=ax, logx=True, logy=True, bins=100, contours=[0.3, 0.5, 0.8], cmap=cmap
    )
    labels.append(setting)
    proxies.append(proxy_line(color=plt.get_cmap(cmap)(0.5)))

ax.legend(proxies, 
          labels,
          facecolor=None,
          frameon=False, 
          bbox_to_anchor=(1,0), 
          loc='lower right', 
          fontsize=12);<jupyter_output><empty_output><jupyter_text>However, when we add other tectonic settings to the diagram, we find something reminiscent of Li's diagram shown above, with a large degree of overlap:<jupyter_code>for setting, ID, cmap in [("Back Arcs", "BAB", "Purples")]:
    df.loc[(df.Class == ID), ["Nb/Yb", "Th/Yb"]].pyroplot.density(
        ax=ax, logx=True, logy=True, bins=100, contours=[0.3, 0.5, 0.8], cmap=cmap
    )
    labels.append(setting)
    proxies.append(proxy_line(color=plt.get_cmap(cmap)(0.5)))
    
ax.legend(proxies, 
          labels,
          facecolor=None,
          frameon=False, 
          bbox_to_anchor=(1,0), 
          loc='lower right', 
          fontsize=12)
ax.figure<jupyter_output><empty_output><jupyter_text>Building and Testing ModelsBelow we'll build a machine tectonic discrimination 'pipeline', including necessary data pre-processing, training and selecting a classifier, and assessing the overall performance of this classification model.There are a variety of types of classifiers we could use to achieve this, and here we've just picked one to work with which is relatively robust and can handle the dimensionality - a support vector classifier (SVC) with a radial basis function. We've included an illustrative example below for a dataset with two different groups of samples (two 'classes'), highlighting some of the differences between some commonly used classifiers.|   ||:--:|| Examples of decision surfaces of different classifier models for the same two-class dataset.| Pre-ProcessingTo make the most of this dataset and build a more precise classifier, we need to make the most of the dataset dimensionality, as using more dimensions gives us more discriminating power. Making a machine learning classifier with these same two dimensions won't get us any further!Most statistical and machine learning algorithms perform best on scaled and approximately normally-distributed data. As we're using major elements to build our classification model, we should **consider using compositional log-ratio transforms** ([Aitchison, 1984]) before we train our model, such that spurious correlation due to closure is minimized (note that we're *not* transforming the isotope ratios here). Also, we should **include a scaling step** to transform our input values such that they fall within $[-1, 1]$, such that magnitude differences between variables are accounted for.[Aitchison, 1984]: https://doi.org/10.1007/BF01029316 "<NAME>., 1984. The statistical analysis of geochemical compositions. Journal of the International Association for Mathematical Geology 16, 531–564. doi: 10.1007/BF01029316" Training a ModelWithin our pipeline, we've specified a range of parameters for the classifier model ('hyperparameters') which it iterates through over a grid to find the most performant model. To ensure that our classifier is robust, we set aside independent subsets of our data for training and testing (used to fit and asses the model, respectively - typically about 30% for testing). We also use cross validation and ensure that separate subsets of the training data are used for fitting and validating models during model selection and tuning processes. Specifically here we used 5-fold stratified cross-validation, which splits the data into five separate representative subsets, and iterates through using each as the validation set. After this process has completed, the held-back test set can be used to independently verify the performance of the classifier model(s). If working as expected and with the desired accuracy and precision, we can extract the best classifier and put it to work on a new dataset! Below we train a classifier using our dataset, print out a quick diagrammatic summary (click for more details) and dig into the overall performance.<jupyter_code>import numpy as np
import pandas as pd
import joblib
import matplotlib.pyplot as plt
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn import set_config
from sklearn.pipeline import FeatureUnion, make_pipeline
set_config(display='diagram')

from pyrolite.util.skl.select import CompositionalSelector
from pyrolite.util.skl.vis import plot_confusion_matrix
from pyrolite.util.skl.pipeline import SVC_pipeline
from pyrolite.util.skl.transform import CLRTransform, LogTransform 
clf_path = "./data/basalt_classifier.pkl"

X = df[[c for c in df.columns if c!='Class']]
y = df['Class']

XX_train, XX_test, yy_train, yy_test = train_test_split(
    X, y, stratify=y, test_size=0.3
)

try:
    gs = joblib.load(clf_path)
except:
    gs = SVC_pipeline(
        verbose=0,
        probability=True,
        scaler=StandardScaler(),
        transform = FeatureUnion([
            ('abundances --> CLR', make_pipeline(CompositionalSelector(), CLRTransform())), 
            ('isotopes', make_pipeline(CompositionalSelector(inverse=True), LogTransform()))
                            ]),
        param_grid=dict(
            svc__gamma=np.logspace(-2.3, -0.85, 5),
            svc__C=np.logspace(-0.5, 2, 5),  # limit C parameter to avoid overfitting
        ),
        n_jobs=1,
        cv=5
    )
    gs.fit(XX_train, yy_train)
    joblib.dump(gs, clf_path)
gs<jupyter_output><empty_output><jupyter_text>Now that we've fit our classifier pipeline, we can extract the best estimator and see how it performed on the test data. While we could simply assess the overall accuracy by checking how many of the test samples it correctly classifies, we can go further and investigate how it does for each of the different tectonic settings within the test dataset. Below we've included a **confusion matrix**, which illustrates the proportion of the samples from a 'true' tectonic setting which are correctly classified (along the diagonal), and also what those which were incorrectly classified were misclassified as (off diagonals along the rows). For example, the back arc basin basalts tend to exhibit poor performance - and we can see that where they are misclassified they're commonly assigned as island arcs, intra-oceanic arcs, continental arcs and MORBs. This stems from geochemical overlap even in 20+ dimensions (as we've seen above in the extreme 2D case), which itself relates to the tectonic proximity and similarity of BAB compared to MORB and island arcs. Similarly, in some cases oceanic plateaus can be difficult to distinguish from continental flood basalts.<jupyter_code>plot_confusion_matrix(gs.best_estimator_, XX_test, yy_test, normalize=True);<jupyter_output><empty_output><jupyter_text>While we've investigated the performance of this classifier a little based on the numbers, let's quickly check what this classifier looks like relative to our Pearce plot. Below we've included a plot colored by the *predicted* tectonic setting for each sample in our dataset, and can see that by exploiting the extra dimensionality, we've managed to untangle some of the overlaps in this space.<jupyter_code>from pyrolite.plot.color import process_color
fig, ax = plt.subplots(1)
pearceThNbYb(ax)
colors = process_color(gs.best_estimator_.predict(X), cmap='tab10', alpha=0.8)['c']
for ID in df.Class.unique():
    subdf = df.loc[df.Class==ID, :]
    subdf[['Nb/Yb', 'Th/Yb']].pyroplot.scatter(c=colors[df.Class==ID], ax=ax, s=8, alpha=0.8)
    
ax.legend([proxy_line(c=coll.get_facecolors()[0], marker='o', ls='none') for coll in ax.collections],
          df.Class.unique(),
          frameon=False,
          facecolor=None)
plt.show()<jupyter_output><empty_output><jupyter_text>There's an extra dimension to this model which we haven't yet explored: uncertainty. While it's relatively straightforward to envisage uncertainty for a predicted variable when it's numeric, here we have a multi-class classification problem. When we built this classifier, we did so using a method which allows for estimation of the relative probability that each sample belongs to each tectonic setting (through 'Platt scaling' of the decision surface; [Platt, 2000]; [Lin, 2007]). While this isn't the most robust way to estimate classification probability, it does provide some useful information about relative certainty.Below we've shown the estimated classification probabilities for the first few samples in our dataset across each of the eight tectonic settings, along with the true tectonic setting. Most samples are classified correctly with high certainties, for others the classification probability is distributed over multiple (typically geologically similar) settings. Where samples are incorrectly classified, it's typically with low certainty.[Platt, 2000]: "" "<NAME>., 2000. Probabilistic Outputs for Support Vector Machines and Comparisons to Regularized Likelihood Methods, in: Advances in Large Margin Classifiers. MIT Press, pp. 61–74."[Lin, 2007]:https://doi.org/10/ckf84q "<NAME>., <NAME>., <NAME>., 2007. A note on Platt’s probabilistic outputs for support vector machines. Mach Learn 68, 267–276."<jupyter_code>samples = np.random.choice(np.arange(XX_test.shape[0]), 8, replace=False)
probs = (
    pd.DataFrame(
        data=gs.best_estimator_.predict_proba(XX_test.iloc[samples]),
        columns=gs.best_estimator_.classes_,
    )
    * 100
)
probs["True Tectonic Setting"] = yy_test.iloc[samples].tolist()
table = probs[probs.columns[::-1]].round(2)
table.style.background_gradient(cmap="Blues", axis=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07976860205922938, 'fraction_numerical': 0.02600573187559707, 'mean_word_length': 4.7395674687785565, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 55, 'https://': 10, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '614609', 'license': 'MIT', 'n_tokens_mistral': 5546, 'n_tokens_neox': 4840, 'n_words': 2234, 'path': '011_classification.ipynb', 'repo_name': 'morganjwilliams/gs2020-diggingdeeper', 'size': '18959'}