01_training_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This is for creating the training data with different numbers of adjacent wells\n",
    "It will write 22 `csv` files to the directory for the varying numbers of adjacent wells. Similar to `01_training_data_generation.ipynb` but with more adjacent wells in the feature set.\n",
    "\n",
    "These 300 well dataset can be downloaded from https://osf.io/a6cwh/ inside the `Training Datasets` folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial.distance import pdist, squareform\n",
    "from sklearn.preprocessing import FunctionTransformer\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Below we define functions that build the dataset. Functions are descriptively named and have full documentation for each function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flatten(container):\n",
    "    \"Flattens lists\"\n",
    "    for i in container:\n",
    "        if isinstance(i, (list, tuple)):\n",
    "            for j in flatten(i):\n",
    "                yield j\n",
    "        else:\n",
    "            yield i\n",
    "\n",
    "def feature_list(no_of_neighbors):\n",
    "    \"\"\"\n",
    "    Creates a list of features given number of adjacent wells\n",
    "    param no_of_neighbors: number of adjacent wells used in feature engineering\n",
    "    \"\"\"\n",
    "    print(\"Getting the features\")\n",
    "    initial = [\"thickness\", \"thickness natural log\", \"thickness power\"]\n",
    "    features = []\n",
    "    for item in initial:\n",
    "        features.append(item)\n",
    "        for i in range(1, no_of_neighbors + 1):\n",
    "            features.append(item + \" neighbor \" + str(i))\n",
    "    features.append([\"x location\", \"y location\", \"class\"])\n",
    "    return list(flatten(features))\n",
    "\n",
    "def data_normalization(dataframe):\n",
    "    \"\"\"\n",
    "    Normalizes the generated stratigraphic data\n",
    "    param dataframe: a pandas dataframe of thicknesses that has been through feature engineering\n",
    "    param location: a pandas dataframe with locations \n",
    "    \"\"\"\n",
    "    print(\"normalizing\")\n",
    "    normalized_df = (dataframe - dataframe.min()) / (dataframe.max() - dataframe.min()).replace(0, 0.00001)\n",
    "    return normalized_df\n",
    "    \n",
    "def rotation(dataframe, j):\n",
    "    \"\"\"\n",
    "    Creates spatial samples and rotates them in the xy plane\n",
    "    param dataframe: dataframe output from stratigraphy generation\n",
    "    param j: controls the rotation of the dataset 0 is no rotation\n",
    "    \"\"\"\n",
    "    x = np.arange(0, 40, 0.1)\n",
    "    y = np.random.randint(0, 10, len(x))\n",
    "    # this is the rotation of the generated data\n",
    "    if j % 0.2 > 0.1:\n",
    "        dataframe[\"ex\"] = x * np.cos(-j / 2) - y * np.sin(-j / 2)\n",
    "        dataframe[\"ey\"] = y * np.cos(-j / 2) - x * np.sin(-j / 2)\n",
    "    else:\n",
    "        dataframe[\"ex\"] = x * np.cos(j / 2) - y * np.sin(j / 2)\n",
    "        dataframe[\"ey\"] = y * np.cos(j / 2) - x * np.sin(j / 2)\n",
    "    return dataframe\n",
    "\n",
    "\n",
    "def missing(dataframe, number_of_layers):\n",
    "    \"\"\"\n",
    "    Inserts missing data into the dataset at random\n",
    "    param dataframe: the dataframe output from the rotation function\n",
    "    param number_of_layers: number of layers to evaluate\n",
    "    \"\"\"\n",
    "    for k in range(200):\n",
    "        dataframe.iloc[\n",
    "            np.random.randint(0, 399), np.random.randint(0, number_of_layers - 1),\n",
    "        ] = 0\n",
    "    return dataframe\n",
    "\n",
    "        \n",
    "def adjacent_wells(dataframe, no_of_neighbors):\n",
    "    \"\"\"\n",
    "    Calculates the adjacent wells and builds the initial dataframe\n",
    "    param dataframe: dataframe output from the missing function\n",
    "    param no_of_neighbors: number of adjacent wells used in feature engineering\n",
    "    \"\"\"\n",
    "    hood = squareform(pdist(dataframe.iloc[:, -2:]))\n",
    "    neighbors = []\n",
    "    for i in enumerate(hood.argsort()[0:, 1 : no_of_neighbors + 1]):\n",
    "        selected = (\n",
    "            dataframe.iloc[hood.argsort()[i[0], 1 : no_of_neighbors + 1], 0:-2]\n",
    "            .stack()\n",
    "            .to_frame()\n",
    "            .T\n",
    "        )\n",
    "        selected.columns = selected.columns.droplevel()\n",
    "        neighbors.append(selected)\n",
    "    frame = pd.concat(neighbors, sort=False)\n",
    "    frame.index = range(len(frame))\n",
    "    neighborhood = pd.concat([dataframe.iloc[:, :-2], frame], axis=1)\n",
    "    return neighborhood\n",
    "    \n",
    "def depth_to_thickness(neighborhood, dataframe):\n",
    "    \"\"\"\n",
    "    Converts the depth dataframe from the adjacent wells function to thicknesses\n",
    "    param neighborhood: dataframe output from `adjacent_wells`\n",
    "    param dataframe: dataframe output from function `missing`\n",
    "    \"\"\"\n",
    "    locations = pd.DataFrame()\n",
    "    df = pd.DataFrame()\n",
    "    thicknesses = neighborhood.diff(axis=1)\n",
    "    thicknesses[thicknesses < 0] = 0\n",
    "    thicknesses.drop(columns=\"zero\", inplace=True)\n",
    "    locations = pd.concat((locations, dataframe.iloc[:, -2:]))\n",
    "    df = pd.concat((df, thicknesses))\n",
    "    return df, locations\n",
    "  \n",
    "    \n",
    "def truncation(smallest, largest, step, names, number_of_layers, j):\n",
    "    \"\"\"\n",
    "    Creates truncated stratal geometries using a min, max, step, names and numbers of layers\n",
    "    param smallest: the smallest integer value for stratigraphy\n",
    "    param largest: the largest integer value for stratigraphy\n",
    "    param step: the size of the step from smallest to largest\n",
    "    param names: names of the layers as strings in a list\n",
    "    param number_of_layers: number of layers to evaluate\n",
    "    param j: float value \n",
    "    \"\"\"\n",
    "    rolling = pd.DataFrame()\n",
    "    j = np.round(j, decimals=3) + 0.5\n",
    "    elevation_random = sorted(\n",
    "        np.random.uniform(smallest, largest, number_of_layers - 1)\n",
    "    )\n",
    "    for i in range(len(names[0 : number_of_layers - 1])):\n",
    "\n",
    "        basement = (\n",
    "            0.001\n",
    "            + (10) * np.sin(1 - np.arange(0, 40, 0.1) / (j * 2) + 0.001)\n",
    "            + np.random.rand(400) / 5\n",
    "        )\n",
    "        elevation = (\n",
    "            np.full(\n",
    "                400,\n",
    "                basement.max()\n",
    "                + np.random.uniform(basement.min() / 2, basement.max() / 64, 1),\n",
    "            )\n",
    "            + np.random.rand(400) / 5\n",
    "        )\n",
    "        topbasement = np.where(basement > elevation, elevation, basement)\n",
    "\n",
    "        rolling[\"zero\"] = topbasement\n",
    "        layer_elevation = (\n",
    "            0.001\n",
    "            + (10) * np.sin(1 - np.arange(0, 40, 0.1) / (j * 2) + 0.001)\n",
    "            + abs(elevation_random[i])\n",
    "            + np.random.rand(400) / 5\n",
    "        )\n",
    "        layer_elevation = np.where(\n",
    "            layer_elevation < basement, basement, layer_elevation\n",
    "        )\n",
    "        layer_elevation = np.where(\n",
    "            layer_elevation > elevation, elevation, layer_elevation\n",
    "        )\n",
    "        rolling[names[i]] = layer_elevation\n",
    "    return rolling\n",
    "\n",
    "def onlap(smallest, largest, step, names, number_of_layers, j):\n",
    "    \"\"\"\n",
    "    Creates onlap stratal geometries using a min, max, step, names and numbers of layers\n",
    "    param smallest: the smallest integer value for stratigraphy\n",
    "    param largest: the largest integer value for stratigraphy\n",
    "    param step: the size of the step from smallest to largest\n",
    "    param names: names of the layers as strings in a list\n",
    "    param number_of_layers: number of layers to evaluate\n",
    "    param j: float value \n",
    "    \"\"\"\n",
    "    rolling = pd.DataFrame()\n",
    "    j = np.round(j, decimals=3) + 0.5\n",
    "    elevation_random = sorted(\n",
    "        np.random.uniform(smallest, largest, number_of_layers - 1)\n",
    "    )\n",
    "    for i in range(len(names[0 : number_of_layers - 1])):\n",
    "        basement = (\n",
    "            0.001\n",
    "            + (10) * np.sin(1 - np.arange(0, 40, 0.1) / (j * 2) + 0.001)\n",
    "            + np.random.rand(400) / 5\n",
    "        )\n",
    "        elevation = (\n",
    "            np.full(\n",
    "                400,\n",
    "                basement.max()\n",
    "                + np.random.uniform(basement.min() / 2, basement.max() / 64, 1),\n",
    "            )\n",
    "            + np.random.rand(400) / 5\n",
    "        )\n",
    "        topbasement = np.where(basement > elevation, elevation, basement)\n",
    "        rolling[\"zero\"] = topbasement\n",
    "        strat_elevation = (\n",
    "            np.full(400, elevation_random[i]) + np.random.rand(400) / 5\n",
    "        )\n",
    "        onlap = np.where(strat_elevation > basement, strat_elevation, basement)\n",
    "        layer_elevation = np.where(onlap > elevation, elevation, onlap)\n",
    "        rolling[names[i]] = layer_elevation\n",
    "    return rolling\n",
    "\n",
    "def horizontal(smallest, largest, step, names, number_of_layers, j):\n",
    "    \"\"\"\n",
    "    Creates onlap stratal geometries using a min, max, step, names and numbers of layers\n",
    "    param smallest: the smallest integer value for stratigraphy\n",
    "    param largest: the largest integer value for stratigraphy\n",
    "    param step: the size of the step from smallest to largest\n",
    "    param names: names of the layers as strings in a list\n",
    "    param number_of_layers: number of layers to evaluate\n",
    "    param j: float value \n",
    "    \"\"\"\n",
    "    rolling = pd.DataFrame()\n",
    "    j = j + 0.5\n",
    "    elevation_random = sorted(\n",
    "        np.random.uniform(smallest, largest, number_of_layers - 1)\n",
    "    )\n",
    "    for i in range(len(names[0 : number_of_layers - 1])):\n",
    "        strat_elevation = (\n",
    "            np.full(400, elevation_random[i]) + np.random.rand(400) / 5\n",
    "        )\n",
    "        basement = strat_elevation - abs(\n",
    "            np.random.uniform(smallest, largest, number_of_layers - 1)\n",
    "            + np.random.rand(400) / 5\n",
    "        )\n",
    "        elevation = (\n",
    "            np.full(400, strat_elevation + elevation_random[i])\n",
    "            + np.random.rand(400) / 5\n",
    "        )\n",
    "        topbasement = np.where(basement > elevation, elevation, basement)\n",
    "        layer_elevation = np.where(\n",
    "            strat_elevation > elevation, elevation, strat_elevation\n",
    "        )\n",
    "        rolling[\"zero\"] = topbasement\n",
    "        rolling[names[i]] = layer_elevation\n",
    "    return rolling\n",
    "\n",
    "def build_feature_engineered_dataset(thickness_df, locations_df):\n",
    "    \"\"\"\n",
    "    Takes the generated thickness dataset and runs feature engineering\n",
    "    param thickness_df: the generated thickness dataset\n",
    "    param locations_df: the generated locations dataset    \n",
    "    \"\"\"\n",
    "    log_transform = FunctionTransformer(np.log, validate=False)\n",
    "    power_transform = FunctionTransformer(lambda x: x ** 10, validate=False)\n",
    "    \n",
    "    logged = pd.DataFrame(log_transform.transform(thickness_df))\n",
    "    powered = pd.DataFrame(power_transform.transform(thickness_df)) \n",
    "    feature_dataset = (\n",
    "        pd.concat([thickness_df, logged, powered, locations_df], axis=1)\n",
    "        .dropna()\n",
    "        .replace(-np.inf, 0)\n",
    "    )\n",
    "    return feature_dataset\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define some global variables to build the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "NAMES = [\"one\", \"two\", \"three\"]  # this creates dummy NAMES for the formations\n",
    "NUMBER_OF_LAYERS = 2  # this is the number of tops you want in your training data\n",
    "SMALLEST = -12\n",
    "LARGEST = 12\n",
    "STEP = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last we run through the different number of neighbors to generate for each dataset and save them to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": "Getting the features\nSTARTING with 399 neighbor\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nnormalizing\nsaving the training data for 399\nDone with 399 neighbors\n"
    }
   ],
   "source": [
    "no_of_neighbors = 399\n",
    "np.random.seed(19)\n",
    "df = pd.DataFrame()\n",
    "df_onlap = pd.DataFrame()\n",
    "locations = pd.DataFrame()\n",
    "elevation_random = sorted(\n",
    "    np.random.uniform(SMALLEST, LARGEST, NUMBER_OF_LAYERS - 1)\n",
    ")\n",
    "flat_features = feature_list(no_of_neighbors)\n",
    "print(f\"STARTING with {no_of_neighbors} neighbor\")\n",
    "trunc_master = pd.DataFrame()\n",
    "onlap_master = pd.DataFrame()\n",
    "horiz_master = pd.DataFrame()\n",
    "\n",
    "for j in np.arange(SMALLEST, LARGEST, STEP):\n",
    "    # first create the truncated dataset\n",
    "    trunc_data = truncation(SMALLEST, LARGEST, STEP, NAMES, NUMBER_OF_LAYERS, j)\n",
    "    trunc_rotated = rotation(trunc_data, j)\n",
    "    trunc_missing = missing(trunc_rotated, NUMBER_OF_LAYERS)\n",
    "    trunc_depth = adjacent_wells(trunc_missing, no_of_neighbors)\n",
    "    trunc_thickness, trunc_locations = depth_to_thickness(trunc_depth, trunc_missing)\n",
    "\n",
    "    # next make the onlap dataset\n",
    "    onlap_data = onlap(SMALLEST, LARGEST, STEP, NAMES, NUMBER_OF_LAYERS, j)\n",
    "    onlap_rotated = rotation(onlap_data, j)\n",
    "    onlap_missing = missing(onlap_rotated, NUMBER_OF_LAYERS)\n",
    "    onlap_depth = adjacent_wells(onlap_missing, no_of_neighbors)\n",
    "    onlap_thickness, onlap_locations = depth_to_thickness(onlap_depth, onlap_missing)\n",
    "\n",
    "    # last make the horizontal dataset\n",
    "    horiz_data = horizontal(SMALLEST, LARGEST, STEP, NAMES, NUMBER_OF_LAYERS, j)\n",
    "    horiz_rotated = rotation(horiz_data, j)\n",
    "    horiz_missing = missing(horiz_rotated, NUMBER_OF_LAYERS)\n",
    "    horiz_depth = adjacent_wells(horiz_missing, no_of_neighbors)\n",
    "    horiz_thickness, horiz_locations = depth_to_thickness(horiz_depth, horiz_missing)\n",
    "\n",
    "    # Now let's run the feature engineering\n",
    "    trunc_non_normal = build_feature_engineered_dataset(trunc_thickness, trunc_locations)\n",
    "    onlap_non_normal = build_feature_engineered_dataset(onlap_thickness, onlap_locations)\n",
    "    horiz_non_normal = build_feature_engineered_dataset(horiz_thickness, horiz_locations)\n",
    "\n",
    "    # Finally we need to normalize the data with our own MinMaxScaler\n",
    "    trunc_normal = data_normalization(trunc_non_normal)\n",
    "    onlap_normal = data_normalization(onlap_non_normal)\n",
    "    horiz_normal = data_normalization(horiz_non_normal)\n",
    "\n",
    "    #Then add class labels\n",
    "    trunc_normal[\"class\"] = 'truncation'  # truncation\n",
    "    onlap_normal[\"class\"] = 'onlap'  # onlap\n",
    "    horiz_normal[\"class\"] = 'horizontal'  # horizontal\n",
    "\n",
    "    trunc_master = pd.concat((trunc_master, trunc_normal))\n",
    "    onlap_master = pd.concat((onlap_master, onlap_normal))\n",
    "    horiz_master = pd.concat((horiz_master, horiz_normal))\n",
    "\n",
    "dataset = pd.concat((trunc_master, onlap_master, horiz_master))\n",
    "dataset.columns = flat_features\n",
    "print(f\"saving the training data for {no_of_neighbors}\")\n",
    "dataset.to_csv(\"stratigraphic_geometry_dataset.csv\")\n",
    "print(f\"Done with {no_of_neighbors} neighbors\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3-final"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}