<jupyter_start><jupyter_text>Cross-validation and hyperparameter tuning for neural networksIn order to do a meaninguful comparison between the decision trees based on the iMeta algorithm, and the neural networ apporach taken in the first paper, we need to run cross validation. Training is comparitively expensive, so we will do fewer splits and fewer param combos, until such time as more parallel processing can be implemented.  Set up toolsStart by importing libraries and specifying the parameters for the experiment.<jupyter_code>import warnings
warnings.filterwarnings('ignore')
import os
import sys
import pathlib
import ipywidgets
import time
import pandas
import numpy
import matplotlib
import matplotlib.pyplot
import sklearn
import sklearn.model_selection
import sklearn.preprocessing
import sklearn.neural_network
import sklearn.metrics
root_repo_dir = pathlib.Path().absolute().parent
sys.path = [os.path.join(root_repo_dir)] + sys.path
import xbt.dataset
from xbt.dataset import XbtDataset, UNKNOWN_STR, cat_output_formatter, check_value_found
from xbt.imeta import imeta_classification, XBT_MAX_DEPTH
from imeta import imeta_classification, XBT_MAX_DEPTH
# Set up some site specific parameters for the notebook
try:
    environment = os.environ['XBT_ENV_NAME']
except KeyError:
    environment = 'pangeo'
root_data_dirs = {
    'MO_scitools': '/data/users/shaddad/xbt-data/',
    'pangeo': '/data/misc/xbt-data/',
}
env_date_ranges = {
    'MO_scitools': (1966,2015),
    'pangeo': (1966,2015)
}<jupyter_output><empty_output><jupyter_text>Define experiment parametersThese are some key hyperparameters we're using for this experiment, such as input features and number of splits.<jupyter_code>input_feature_names = ['country','max_depth', 'year']
target_feature = 'instrument'
splitter_features = ['year', 'instrument']
unseen_feature = 'cruise_number'
num_unseen_splits = 3
unseen_fraction = 1.0 / num_unseen_splits
num_train_splits = 3
train_fraction = 1.0 - (1.0 / num_train_splits)
unseen_fold_name = 'unseen_fold'
experiment_name = 'nb_mlp_cv_country'
classifier_class = sklearn.neural_network.MLPClassifier
classifier_opts = {'max_iter': 1000,
                  'hidden_layer_sizes':[10, 10,],
                   'alpha': 1e-4,
                  }
classifier_name = 'MLP'
suffix='country'
classifier_param_grid = {
    'hidden_layer_sizes':[(5, 5,), (10, 10,)]
},
cv_metrics = ['recall_micro', 'precision_micro', 'f1_micro']
hp_tune_metric = 'recall_micro'
# Set up some dataset specific parameters
root_data_dir = root_data_dirs[environment]
year_range = env_date_ranges[environment]
input_dir_name = 'csv_with_imeta'
exp_out_dir_name = 'experiment_outputs'
xbt_input_dir = os.path.join(root_data_dir, input_dir_name)
xbt_output_dir = os.path.join(root_data_dir, exp_out_dir_name, experiment_name)
# create the output for this experiment if it doesn't exist
if not os.path.isdir(xbt_output_dir):
    os.makedirs(xbt_output_dir)
output_fname_template = 'xbt_output_{classifier}_{suffix}.csv'
result_fname_template = 'xbt_metrics_{classifier}_{suffix}.csv'<jupyter_output><empty_output><jupyter_text>Running the cross validation experimentThere are multiple levels to the experiment we want to run. The first to train on part of the data and then test on another part. We want to divide up the profiles by cruise, so the test data consists of cruises that were not seen in the training data. We then want to train a model for each of the possible train/test splits. We divide the model into "folds", and then each splits uses data in fold n as the test set and ther rest as the train, as is standard in machine learning. In this case calculate the fold labels manually rather than using standard functions so we can divide based on cruise number and assign profiles to a fold based on which cruise they come from.For the training on eeach split, we also then want to use the best hyperparameters. So for each split we run a hyperparameter tuning scheme, ion this case a grid search through the request values of the hyper parameters. The hyper parameters used here are:* max depth (of the decision tree, not the temperature probe in this case!!)* min_samples_leaf -  a parameter controlling how frequently a leaf node splits. Lower numbers redces over-fitting tendencies* citerion - different criteria for scoring predictions while training.The bottom level is that for each split and set of hyperparameters, we also want to run cross valdation within that training set, to check we have robust results. So we haveter and inner cross validation run together with hper parameter tuning to ensure we have robust results that are as good as we can get, and that demonstrate how well the trainined algorithm will generalise to the unlabelled data. More information on the approach taken here and the spceific classes and function used can found in the following links:scikit-learn API documentation:* https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.htmlsklearn.tree.DecisionTreeClassifier* https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_validate.html* https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GroupKFold.html* https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.htmldiscussion of cross-validationand hyperparameters for decision trees:* https://scikit-learn.org/stable/modules/cross_validation.html* https://stackoverflow.com/questions/60996995/use-groupkfold-in-nested-cross-validation-using-sklearn* https://scikit-learn.org/stable/auto_examples/model_selection/plot_nested_cross_validation_iris.html* https://scikit-learn.org/stable/modules/tree.html<jupyter_code>%%time
xbt_full_dataset = XbtDataset(xbt_input_dir, year_range)
%%time
xbt_labelled = xbt_full_dataset.filter_obs({'labelled': 'labelled'})
_ = xbt_labelled.get_ml_dataset(return_data=False)
_ = xbt_labelled.filter_features(['instrument','model','manufacturer']).encode_target(return_data = False)
%%time
xbt_labelled.generate_folds_by_feature('cruise_number', num_unseen_splits, unseen_fold_name)
X_labelled = xbt_labelled.filter_features(input_feature_names).get_ml_dataset()[0]
y_labelled = xbt_labelled.filter_features([target_feature]).get_ml_dataset()[0]

group_cv1 = sklearn.model_selection.GroupKFold(n_splits=num_unseen_splits)

clf_dt1 = classifier_class(**classifier_opts)
clf_gridcv1 = sklearn.model_selection.GridSearchCV(clf_dt1,
                                                   param_grid = classifier_param_grid,                                          
                                                   scoring=hp_tune_metric,
                                                   cv=num_train_splits,
                                                  )
%%time
scores = sklearn.model_selection.cross_validate(
    clf_gridcv1,
    X_labelled, y_labelled, 
    groups=xbt_labelled[unseen_fold_name], 
    cv=group_cv1,
    return_estimator=True,
    return_train_score=True,
    scoring=cv_metrics,
    n_jobs=-1,
)<jupyter_output>CPU times: user 557 ms, sys: 645 ms, total: 1.2 s
Wall time: 1h 33min 25s<jupyter_text>Plotting results of cross-validation and hyperparameter tuningIn this section we will take a look at the results from the cross validation and hyperparameter tuning of the decision tree to the XBT labelled data. To gain greater understanding, we seprately score the train and test sets for each year of the data and each class of instrument label. This is calculated separately for eeach of the estimators trained for the outer cross-validation (based on cruise numbers). We calculate the recall, precision and f1 metrics both for each class and averaged over all classes.<jupyter_code>def generate_metrics(clf, xbt_ds, data_label, target_feature, input_feature_names, year_range):
    metric_list = []
    for year in range(year_range[0],year_range[1]):
        X_year = xbt_ds.filter_obs({'year': year}, ).filter_features(input_feature_names).get_ml_dataset()[0]
        y_year = xbt_ds.filter_obs({'year': year} ).filter_features([target_feature]).get_ml_dataset()[0]

        y_res_year = clf.predict(X_year)
        cats = list(xbt_ds._feature_encoders[target_feature].classes_)
        prec_year, recall_year, f1_year, support_year = sklearn.metrics.precision_recall_fscore_support(
            y_year, y_res_year, average='micro', labels=range(0,len(cats)))
        prec_cat, recall_cat, f1_cat, support_cat = sklearn.metrics.precision_recall_fscore_support(
            y_year, y_res_year, labels=range(0,len(cats)))

        column_template = '{metric}_{data}_{subset}'
        metric_dict = {'year': year,
                       column_template.format(data=data_label, metric='precision', subset='all'): prec_year,
                       column_template.format(data=data_label, metric='recall', subset='all'): recall_year,
                       column_template.format(data=data_label, metric='f1', subset='all'): f1_year,
                      }

        metric_dict.update({column_template.format(data=data_label, metric='precision', subset=cat): val for cat, val in zip(cats, prec_cat)})
        metric_dict.update({column_template.format(data=data_label, metric='recall', subset=cat): val for cat, val in zip(cats, recall_cat)})
        metric_dict.update({column_template.format(data=data_label, metric='f1', subset=cat): val for cat, val in zip(cats, f1_cat)})
        metric_dict.update({column_template.format(data=data_label, metric='support', subset=cat): val for cat, val in zip(cats, support_cat)})
        metric_list += [metric_dict]

    metrics_df = pandas.DataFrame.from_records(metric_list)
    return metrics_df 
%%time
metrics_list = {}
for split_num, estimator in enumerate(scores['estimator']):
    xbt_train = xbt_labelled.filter_obs({unseen_fold_name: split_num}, mode='exclude')
    xbt_test = xbt_labelled.filter_obs({unseen_fold_name: split_num}, mode='include')
    train_metrics = generate_metrics(estimator, xbt_train, 'train_{0}'.format(split_num), target_feature, input_feature_names, year_range )
    test_metrics = generate_metrics(estimator, xbt_test, 'test_{0}'.format(split_num), target_feature, input_feature_names, year_range)
    metrics_list[split_num] = pandas.merge(train_metrics, test_metrics, on='year')

fig_results = matplotlib.pyplot.figure('xbt_results',figsize=(25,15))
for label1, metrics1  in metrics_list.items():
    ax_precision = fig_results.add_subplot(3,5,label1 +1, title='precision split {0}'.format(label1))
    ax_recall = fig_results.add_subplot(3,5,label1 + 1 + 5 * 1, title='recall split {0}'.format(label1))
    ax_f1 = fig_results.add_subplot(3,5,label1 + 1 + 5 * 2, title='f1 split {0}'.format(label1))
    metrics1.plot.line(ax=ax_precision, x='year', y=[f'precision_train_{label1}_all',f'precision_test_{label1}_all'], color=['b', 'r'], ylim=(0.4,1.0))
    metrics1.plot.line(ax=ax_recall, x='year', y=[f'recall_train_{label1}_all',f'recall_test_{label1}_all'], color=['b', 'r'], ylim=(0.4,1.0))
    metrics1.plot.line(ax=ax_f1, x='year', y=[f'f1_train_{label1}_all',f'f1_test_{label1}_all'], color=['b', 'r'], ylim=(0.4,1.0))<jupyter_output><empty_output><jupyter_text>Also of interest is what hyper parameters were chosen by the GridCV hyperparameter tuning process. The following listing shows which parameters result in the best results for each split. We see that limiting the number of nodes resulted in better performance and restricting how nodes could split (min_samples_leaf) also varied based o the split. Varying these parameters results in different classifications so we get some profiles with different class assignments between different classifiers.<jupyter_code>[est1.best_params_ for est1 in scores['estimator']]
metrics_df_merge = None
for label1, metrics1  in metrics_list.items():
    if metrics_df_merge is None:
        metrics_df_merge = metrics1
    else:
        metrics_df_merge = pandas.merge(metrics_df_merge, metrics1)
results_output_path = os.path.join(xbt_output_dir,
                                   result_fname_template.format(classifier=classifier_name, 
                                                                suffix=suffix))
metrics_df_merge.to_csv(results_output_path)<jupyter_output><empty_output><jupyter_text>Outputing the dataTo filter based on what profiles we can use for predicting, we need some way of checking each profile we create these checkers from the labelled dataset, because the subset of data that was used for training determines what subset is valid for prediction. For example, if a country is not present in the training data, then the prediction function won't be able to handle that profile to predict a probe model and manufacturer. Profiles that are not handled by the trained classifier will get the label "UNKNOWN".<jupyter_code># checker functions check each element of the profile metadata that could be a problem. The checkers are constructed from the labelled data subset.
checkers_labelled = {f1: c1 for f1, c1 in xbt_labelled.get_checkers().items() if f1 in input_feature_names}
result_feature_template = 'instrument_res_dt_split{split_num}'
%%time
xbt_predictable = xbt_full_dataset.filter_predictable(checkers_labelled)
%%time
result_feature_names = []
# classifications_df = None
for split_num, estimator in enumerate(scores['estimator']):
    res_name = result_feature_template.format(split_num=split_num)
    result_feature_names += [res_name]
    res_ml1 = estimator.predict(xbt_predictable.filter_features(input_feature_names).get_ml_dataset()[0])
    res2 = list(xbt_labelled._feature_encoders['instrument'].inverse_transform(res_ml1).reshape(-1))
    xbt_predictable.xbt_df[res_name] = res2

def imeta_instrument(row1):
    return 'XBT: {t1[0]} ({t1[1]})'.format(t1=imeta_classification(row1))
%%time
imeta_output_flags = []
for res_name in result_feature_names:
    flag_name = f'imeta_applied_{res_name}'
    imeta_output_flags += [flag_name]
    xbt_predictable.xbt_df[flag_name] = 0
    xbt_predictable.xbt_df.loc[xbt_predictable.xbt_df[xbt_predictable.xbt_df[res_name].isnull()].index, flag_name] = 1
    xbt_predictable.xbt_df[flag_name] = xbt_predictable.xbt_df[flag_name].astype('int8')
    xbt_predictable.xbt_df.loc[xbt_predictable.xbt_df[xbt_predictable.xbt_df[res_name].isnull()].index, res_name] = \
        xbt_predictable.xbt_df[xbt_predictable.xbt_df[res_name].isnull()].apply(imeta_instrument, axis=1)


%%time
fv_dict = {rfn1: UNKNOWN_STR for rfn1 in result_feature_names}
fv_dict.update({f1: 1 for f1 in imeta_output_flags})
xbt_full_dataset.merge_features(xbt_predictable, result_feature_names + imeta_output_flags,
                               fill_values = fv_dict,
                               feature_encoders={rfn1: xbt_labelled._feature_encoders['instrument'] for rfn1 in result_feature_names},
                               target_encoders={rfn1: xbt_labelled._target_encoders['instrument'] for rfn1 in result_feature_names},                                
                                output_formatters={rfn1: cat_output_formatter for rfn1 in result_feature_names})
%%time
xbt_unknown_inputs = xbt_full_dataset.filter_obs({'classification_quality_flag': 0})
imeta_instrument_fallback = xbt_unknown_inputs.xbt_df.apply(imeta_instrument, axis=1)
for res_name in result_feature_names:
    xbt_full_dataset.xbt_df.loc[xbt_unknown_inputs.xbt_df.index, res_name] = imeta_instrument_fallback<jupyter_output>CPU times: user 335 ms, sys: 20.1 ms, total: 355 ms
Wall time: 351 ms<jupyter_text>The classifier is unable to predict a class for some of the profiles. This can be either because the classifier cannot produce a result, or because the country is not present in the labelled dataset. In these cases iMeta is used as a fallback. This plot shows for what proportion of profiles this occurs.<jupyter_code>pandas.DataFrame.from_records([dict(zip(['ML tree','imeta'] ,xbt_full_dataset.xbt_df[imf1].value_counts().values)) for imf1 in imeta_output_flags]).plot(kind='bar',stacked=True,title='Number of profiles falling back on iMeta', figsize=(12,12))<jupyter_output><empty_output><jupyter_text>Of interest in our classifier ensemble is how many profiles are labelled differently by different classifiers. Here we count the different labels for each profile. We say that the ajority get only 1 label, but a substantial minority get 2 ior 3 labels, and some get as many as 5 (a different classification from each trained instance of the classifier!). For those profile which have been assigned more that one label, we can assign probabilities based on the number of classifiers "voting" for each label. This is not a real probability, but give some idea of the confidence we have in assigned labels. It could also be used as a distribution from which to sample for an esemble temperature dataset.<jupyter_code>%%time
count1 = xbt_full_dataset.filter_features(result_feature_names).xbt_df.apply(lambda row1: len(set(row1.values)), axis='columns')
count1.value_counts()
%%time
vote_count = numpy.zeros([xbt_full_dataset.shape[0], len(xbt_full_dataset._feature_encoders[result_feature_names[0]].classes_[0])],dtype=numpy.float64)
for res_name in result_feature_names:
    vote_count += xbt_full_dataset.filter_features([res_name]).get_ml_dataset()[0]
    
vote_count /= float(len(result_feature_names))
%%time
vote_dict = {f'instrument_probability_{cat1}': vote_count[:,ix1] for ix1, cat1 in enumerate(xbt_full_dataset._feature_encoders['instrument'].classes_[0])}
vote_dict.update({'id': xbt_full_dataset['id']})
vote_df = pandas.DataFrame(vote_dict)
xbt_full_dataset.xbt_df = xbt_full_dataset.xbt_df.merge(vote_df, on='id')
%%time
xbt_full_dataset.output_data(os.path.join(xbt_output_dir, 
                                          output_fname_template.format(classifier=classifier_name,
                                                                       suffix=suffix)),
                             target_features=[])<jupyter_output>CPU times: user 1min 29s, sys: 1.65 s, total: 1min 30s
Wall time: 1min 46s<jupyter_text>Summary for ML tree with cross-validation and hyperparameter tuning (under construction)<jupyter_code>import sklearn.ensemble
vote_clf1 = sklearn.ensemble.VotingClassifier([s1.best_estimator_ for s1 in scores['estimator']], voting='soft')
y_res_labelled = vote_clf1.transform(X_labelled)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0897087589654423, 'fraction_numerical': 0.012714627254944577, 'mean_word_length': 4.8896, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 34, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '473859', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5777, 'n_tokens_neox': 5526, 'n_words': 1505, 'path': 'notebooks/first_ml_cross_validate.ipynb', 'repo_name': 'uk-gov-mirror/MetOffice.XBTs_classification', 'size': '18404'}