<jupyter_start><jupyter_text>Load training data<jupyter_code>data = pd.read_csv('train.csv', sep=';')
data.sample(10)<jupyter_output><empty_output><jupyter_text>Inspecting available logs and interpretations The data contains the metadata columns* WELL: well name* DEPTH_MD: measured depth* X_LOC: UTM X coordinate * Y_LOC: UTM Y coordinate* Z_LOC: DEPTH* GROUP: NPD lithostratigraphy group* FORMATION: NPD lithostratgraphy formationThe data contains the well log curves. For example the following:* BS: Bit Size* CALI: Caliper* RDEP: Deep Resistivity* RHOB: Bulk Density* GR: Raw gamma data* SGR: Spectral Gamma Ray* RMED: Medium Resistivity* ROP: Rate of Penetration* NPHI: Neutron Porosity* PEF: Photoelectric Absorption Factor* RSHA: Shallow Resistivity* DTS: Sonic (Sheer Slowness)* DTC: Sonic (Compressional Slowness)See contest page for full list.as well as the interpretation* FORCE_2020_LITHOFACIES_LITHOLOGY: lithology class label* FORCE_2020_LITHOFACIES_CONFIDENCE: confidence in lithology interpretation (1: high, 2: medium, 3: low) There is a total of 83 wells in the training dataset. The `WELL` column is included so the data can be separated per well. This will be necessary for any "windowed" approaches where we use non-local information in the ML features.<jupyter_code>wells = data['WELL'].unique()
wells
len(wells)<jupyter_output><empty_output><jupyter_text>The LITHOLOGY column The LITHOLOGY column contains the interpretation. Each number corresponds to a lithology class, given by the following list:<jupyter_code>lithology_keys = {30000: 'Sandstone',
                 65030: 'Sandstone/Shale',
                 65000: 'Shale',
                 80000: 'Marl',
                 74000: 'Dolomite',
                 70000: 'Limestone',
                 70032: 'Chalk',
                 88000: 'Halite',
                 86000: 'Anhydrite',
                 99000: 'Tuff',
                 90000: 'Coal',
                 93000: 'Basement'}
data['FORCE_2020_LITHOFACIES_LITHOLOGY'].unique()<jupyter_output><empty_output><jupyter_text>Lithologies are not equally likely. This is reflected in this dataset.<jupyter_code>counts = data['FORCE_2020_LITHOFACIES_LITHOLOGY'].value_counts()
names = []
percentage = []
N = data['FORCE_2020_LITHOFACIES_LITHOLOGY'].shape[0]
for item in counts.iteritems():
    names.append(lithology_keys[item[0]])
    percentage.append(float(item[1])/N*100)
fig, ax = plt.subplots(1, 1, figsize=(14, 7))
ax.bar(x=np.arange(len(names)), height=percentage)
ax.set_xticklabels(names, rotation=45)
ax.set_xticks(np.arange(len(names)))
ax.set_ylabel('Lithology presence (\%)')<jupyter_output>d:\anaconda3\envs\force\lib\site-packages\ipykernel_launcher.py:10: UserWarning: FixedFormatter should only be used together with FixedLocator
  # Remove the CWD from sys.path while we load stuff.<jupyter_text>Missing logs **Notice that some of the curves are NaN. This is an imporant aspect of this dataset and of this competition.** The only log columns that are **guaranteed** to be present are WELL, DEPT, and GR. All other logs can (and will) be missing from some parts of the dataset.<jupyter_code>data.isna().any()
data.isna().sum()<jupyter_output><empty_output><jupyter_text>The following plot shows how large a percentage of training data wells contain at least some depth interval with the given logs. As you can see a couple of logs like RT and DTS are present in less than half of the training wells. **Remember, the test data will have a similar distribution of available logs**<jupyter_code>occurences = np.zeros(25)
for well in data['WELL'].unique():
    occurences += data[data['WELL'] == well].isna().all().astype(int).values[2:-2]
fig, ax = plt.subplots(1, 1, figsize=(14, 7))
ax.bar(x=np.arange(occurences.shape[0]), height=(data.WELL.unique().shape[0]-occurences)/data.WELL.unique().shape[0]*100.0)
ax.set_xticklabels(data.columns[2:-2], rotation=45)
ax.set_xticks(np.arange(occurences.shape[0]))
ax.set_ylabel('Well presence (\%)')<jupyter_output>d:\anaconda3\envs\force\lib\site-packages\ipykernel_launcher.py:6: UserWarning: FixedFormatter should only be used together with FixedLocator<jupyter_text>Look at some log curves for a couple of wells<jupyter_code>well_data = data[data['WELL'] == wells[0]]
fig, axs = plt.subplots(1, len(well_data.columns)-9, figsize=(14, 10), sharey=True)
for ic, col in enumerate(set(well_data.columns)-set(['DEPTH_MD', 'FORCE_2020_LITHOFACIES_LITHOLOGY',
                                                     'FORCE_2020_LITHOFACIES_CONFIDENCE', 'WELL', 'GROUP', 'FORMATION',
                                                    'X_LOC', 'Y_LOC', 'Z_LOC'])):
    axs[ic].plot(well_data[col], well_data['DEPTH_MD'])
    axs[ic].set_xlabel(col)
axs[0].set_ylim(well_data['DEPTH_MD'].values[-1], well_data['DEPTH_MD'].values[0])
well_data = data[data['WELL'] == wells[57]]
fig, axs = plt.subplots(1, len(well_data.columns)-9, figsize=(14, 10), sharey=True)
for ic, col in enumerate(set(well_data.columns)-set(['DEPTH_MD', 'FORCE_2020_LITHOFACIES_LITHOLOGY',
                                                     'FORCE_2020_LITHOFACIES_CONFIDENCE', 'WELL', 'GROUP', 'FORMATION',
                                                    'X_LOC', 'Y_LOC', 'Z_LOC'])):
    axs[ic].plot(well_data[col], well_data['DEPTH_MD'])
    axs[ic].set_xlabel(col)
axs[0].set_ylim(well_data['DEPTH_MD'].values[-1], well_data['DEPTH_MD'].values[0])<jupyter_output><empty_output><jupyter_text>Preparing data for training<jupyter_code>X = data.drop(columns=['WELL', 'FORCE_2020_LITHOFACIES_LITHOLOGY', 'FORCE_2020_LITHOFACIES_CONFIDENCE'])<jupyter_output><empty_output><jupyter_text>As a simple approach we replace nan with zero and hope for the best<jupyter_code>X.fillna(0, inplace=True)<jupyter_output><empty_output><jupyter_text>We should convert lithology labels to numbers 0-11 before training<jupyter_code>y = data['FORCE_2020_LITHOFACIES_LITHOLOGY']
lithology_numbers = {30000: 0,
                 65030: 1,
                 65000: 2,
                 80000: 3,
                 74000: 4,
                 70000: 5,
                 70032: 6,
                 88000: 7,
                 86000: 8,
                 99000: 9,
                 90000: 10,
                 93000: 11}
y = y.map(lithology_numbers)<jupyter_output><empty_output><jupyter_text>Lithologies are not equally likely. Probably smart to take this into account by computing class weights<jupyter_code>np.bincount(y)<jupyter_output><empty_output><jupyter_text>Shuffle the data<jupyter_code>X, y = shuffle(X, y, random_state=42)
X<jupyter_output><empty_output><jupyter_text>Select only the columns you want to use<jupyter_code>X.columns
X = X[['DEPTH_MD', 'CALI', 'RSHA', 'RMED', 'RDEP', 'RHOB', 'GR', 'NPHI', 'PEF', 'DTC', 'SP', 'BS']]<jupyter_output><empty_output><jupyter_text>Scale the features<jupyter_code>scaler = StandardScaler()
X.columns
X = scaler.fit_transform(X)<jupyter_output><empty_output><jupyter_text>Split into train and test data<jupyter_code>X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)<jupyter_output><empty_output><jupyter_text>Applying the custom scoring function Apart from the missing log curves, **one of the key challenges of this machine learning competition is that you are being scored using a custom scoring function**. A prediction is scored using the penalty function$$ S = - \frac{1}{N}\sum_{i=0}^N A_{ \hat{y}_i y_i} $$where $N$ is the number of samples, $y_i$ is prediction for sample $i$, $\hat{y}_i$ is the true target for sample $i$, and $A$ is the penalty matrix given by the followingIn short, the value of the matrix A at row i and column j is the **penalty** given by guessing lithology number i when the correct label is lithology numner j. Notice that the diagonal consists of zeros, no penalty is given for correct predictions.<jupyter_code>A = np.load('penalty_matrix.npy')
A<jupyter_output><empty_output><jupyter_text>The scoring function can be implemented as follows. If used as a loss function in a neural network framwork, it is recommended to implement a more efficient function based on tensor products.<jupyter_code>def score(y_true, y_pred):
    S = 0.0
    y_true = y_true.astype(int)
    y_pred = y_pred.astype(int)
    for i in range(0, y_true.shape[0]):
        S -= A[y_true[i], y_pred[i]]
    return S/y_true.shape[0]<jupyter_output><empty_output><jupyter_text>**IMPORTANT: The implementation above assumes zero-indexed y_true and y_pred. Be careful not to mix zero-indexed and 1-indexed class labels.** Let us see what happens when we test the scoring function on some corner cases:  Lazily guessing SHALE for everything gives an OK score due to the imbalance of the dataset<jupyter_code>score(y_test.values, np.ones(y_test.shape)*2)<jupyter_output><empty_output><jupyter_text>Doing the same for an uncommon class, like COAL, not so much..<jupyter_code>score(y_test.values, np.ones(y_test.shape)*11)<jupyter_output><empty_output><jupyter_text>Of course getting it completely right gives the highest possible score, 0<jupyter_code>score(y_test.values, y_test.values)<jupyter_output><empty_output><jupyter_text>Example prediction on open test dataset  Train an oversimplified dummy model for the sake of the exampleIn order to give examples of how to submit predictions and reproducible code, we will train a RandomForestClassifier on the full feature set where we lazily have set all missing values (NaNs) to zero. This is an example only, following along this approach for the competition is not recommended. However, the template below can be used for generating predictions for the weekly leaderboard from your own model, and as an example on how to ecapsulate your model for final scoring where you have to submit code.Apart from the inherent flaws of the random forest classifier, this example approach has a number of issues and simplifications. For example:* Data is not separated per well. Normalization and shuffling is done across the whole dataset.* The penatly matrix is not taken into account, it is treated as a normal classification problem.* Missing logs are treated as zeros* Depth and formation tops not taken into account* Same model is used for when you only have gamma ray logs as when you have the full set of logs<jupyter_code>model = RandomForestClassifier(n_estimators=50, random_state=42, class_weight='balanced')
model.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Test model<jupyter_code>y_pred_test = model.predict(X_test)
score(y_test.values, y_pred_test)<jupyter_output><empty_output><jupyter_text>Save model and scaler to files<jupyter_code>pickle.dump(model, open('model.pkl', 'wb'))
pickle.dump(scaler, open('scaler.pkl', 'wb'))<jupyter_output><empty_output><jupyter_text>In order to make it easy to run inference on new data, such as the open test dataset, it is recommended to wrap the initialization of the model, all preprocessing steps including normalization, treatment of NaNs and removals of unused columns in a class or pipeline. There are many ways of doing this, one example being a class like below:<jupyter_code>class Model(object):
    def __init__(self, model_file, scaler_file):
        # Load pre-trained model from file
        self.model = pickle.load(open(model_file, 'rb'))
        # Load a "pre-trained" scaler from file
        self.scaler = pickle.load(open(scaler_file, 'rb'))
        
    def _preprocess(self, features):
        # Method to be run before inference. Contains things like
        # stripping unwanted columns, replacing NaNs, and scaling 
        # or normalizing data
        features.fillna(0, inplace=True)
        features = features[['DEPTH_MD', 'CALI', 'RSHA', 'RMED', 'RDEP', 'RHOB', 'GR', 'NPHI', 'PEF', 'DTC', 'SP', 'BS']]
        return self.scaler.transform(features)
        
    def predict(self, features):
        # This function should be able to take in features in their
        # raw, unprocessed form as read from the file test.csv and
        # return predictions as an array integers of the same length
        X = self._preprocess(features)
        return self.model.predict(X)<jupyter_output><empty_output><jupyter_text>The model can then be initialized using persisted data in the form of saved model file (or models files), scaler files, and other persisted data.<jupyter_code>model = Model('model.pkl', 'scaler.pkl')<jupyter_output><empty_output><jupyter_text>Example submission: weekly leaderboardFor the weekly leaderboard that will start September 1st, you are asked to submit predictions on the supplied open dataset in the file `test.csv`. This file contains all the columns that are in the training data except LITHOLOGY and LITHOLOGY_CONFIDENCE. Notice that the test data will also have missing columns. It is therefore important that your inference pipeline has the treatment of these missing columns as part of the pre-processing step.<jupyter_code>open_test_features = pd.read_csv('test.csv', sep=';')
open_test_features.head()
open_test_features.columns
test_prediction = model.predict(open_test_features)
test_prediction<jupyter_output><empty_output><jupyter_text>**Important:** Before submitting we must map these predictions back to the lithology labels<jupyter_code>category_to_lithology = {y:x for x,y in lithology_numbers.items()}
test_prediction_for_submission = np.vectorize(category_to_lithology.get)(test_prediction)<jupyter_output><empty_output><jupyter_text>For submission to the weekly leaderboard, write the predictions as a csv file- single column (no index column or similar)- header: 'lithology'- exactly the same number of rows as the raw test features- use the original lithology labels, 65000, 30000, ..., NOT the encoded ones 0, 1, 2, ...Example:```lithology65000650003000030000300003000030000650306503065030...```<jupyter_code>test_prediction_for_submission
np.savetxt('test_predictions.csv', test_prediction_for_submission, header='lithology', comments='', fmt='%i')<jupyter_output><empty_output><jupyter_text>Submitting code for final scoringAt the end of the competition, the highest performing teams are invited to submit code for blind scoring on a hidden test dataset. The winner will be the team with the best score on the hidden test dataset, not the weekly leaderboard. In order to do this, you must submit the following in a zip-file:* Code to do the inference* Any persisted data in the form of pre-trained model or models, scalers, and more* Information about your environmentThere are many ways of doing this. In the example above, the code to do the inference would be the `Model` class in a python-file, the persisted data would be the files `model.pkl` and `scaler.pkl`. Your python environment needs to be replicable by the judges of the competition. There are many ways of doing this. Examples include dumping pip package versions to a file<jupyter_code>!pip freeze > requirements.txt<jupyter_output><empty_output><jupyter_text>If you are using Anaconda, then you can export an envionment file as follows<jupyter_code>!conda env export > environment.yml<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09646217039883398, 'fraction_numerical': 0.022525506823903536, 'mean_word_length': 4.7924021488871835, 'pattern_counts': {'":': 0, '<': 119, '<?xml version=': 0, '>': 121, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '30', 'dir': 'jupyter-structured-clean-dedup', 'id': '189903', 'license': 'MIT', 'n_tokens_mistral': 4999, 'n_tokens_neox': 4624, 'n_words': 1526, 'path': 'starter_notebook.ipynb', 'repo_name': 'nguyen-tr-le/TheFORCE2020-WellLogs-TSClassification', 'size': '15094'}