infinite_pool_qbdc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Query by Monte Carlo Dropout Committee\n",
    "### Import packages\n",
    "**Keras**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.models import Sequential, load_model\n",
    "from keras.layers.core import Dense, Activation, Lambda\n",
    "from keras.layers import Dropout\n",
    "from keras.callbacks import History, ModelCheckpoint, EarlyStopping, ReduceLROnPlateau\n",
    "from keras.utils import to_categorical\n",
    "from keras import regularizers\n",
    "from keras import initializers\n",
    "import keras.backend as K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data Science**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "import sklearn.metrics as metrics\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score, log_loss\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Other**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "from math import floor\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Functions\n",
    "**Plot Projection**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_projection(data, classification, x=0, y=1, bins=100):\n",
    "    # Calculate ranges and bins\n",
    "    xmin, xmax = np.amin(data[:,x]), np.amax(data[:,x])\n",
    "    ymin, ymax = np.amin(data[:,y]), np.amax(data[:,y])\n",
    "    xbins = np.linspace(xmin, xmax, bins)\n",
    "    ybins = np.linspace(ymin, ymax, bins)\n",
    "    # Calculate two histograms\n",
    "    allowed, _, _ = np.histogram2d(data[classification==1.0, x], data[classification==1.0, y], [xbins, ybins])\n",
    "    excluded, _, _ = np.histogram2d(data[classification==0.0, x], data[classification==0.0, y], [xbins, ybins])\n",
    "    # Calculate map\n",
    "    mapping = allowed / (allowed + excluded)\n",
    "    mapping = mapping.T\n",
    "    mapping = np.flipud(mapping)\n",
    "    # Plot\n",
    "    f, a = plt.subplots(1,1,figsize=(8,8))\n",
    "    a.matshow(mapping, extent=(xmin, xmax, ymin, ymax), cmap=\"seismic_r\")\n",
    "    print(\"({}, {})\".format(np.amin(mapping),np.amax(mapping)))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Normalisation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mins = np.array([-3999.99, -4000.0, 200.0, 90.879999, 91.419997, 90.02, 90.250001, 200.02999, 102.65, 200.00984, 264.01999, 200.0199, 100.06994, -7998.919, -3999.97, -3999.97, -4381.715, 1364.5547, 2.6724445])\n",
    "maxs = np.array([3999.95, 3999.83, 3999.98, 4000.0, 3999.98, 3999.99, 4000.0, 3999.9899, 3999.9999, 3999.98, 3999.9999, 3999.97, 4000.0, 7993.83, 3999.95, 3999.94, 4132.4793, 37569812.0, 66.371989])\n",
    "\n",
    "def normalize(data):\n",
    "    # Normalize data\n",
    "    mu = (maxs + mins)/1.5\n",
    "    sigma = (maxs - mins)/np.sqrt(12)\n",
    "    data = (data - mu)/sigma\n",
    "    # Return\n",
    "    return data\n",
    "\n",
    "def undo_normalize(data):\n",
    "    # Normalize data\n",
    "    mu = (maxs + mins)/2.0\n",
    "    sigma = (maxs - mins)/np.sqrt(12)\n",
    "    data = data*sigma + mu\n",
    "    # Return\n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Oracle**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oracle_model = load_model(\"/archive/projects/activelearning/susyai_nn/fifth/model.hdf5\")\n",
    "def oracle(data):\n",
    "    return 1.0*(oracle_model.predict(data) > 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Generate data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_data(N, normalized=True):\n",
    "    X = np.random.rand(N,19)\n",
    "    X *= (maxs - mins)\n",
    "    X += mins\n",
    "    Xnormed = normalize(X)\n",
    "    y = oracle(Xnormed).astype(np.float)\n",
    "    if normalize:\n",
    "        return (Xnormed,y)\n",
    "    return (X,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Train Model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model(Xtrain, ytrain, Xval, yval, epochs=1000, verbose=0):\n",
    "    model_name = \"qmcdc_susyai.hdf5\"\n",
    "    # Architecture configuration\n",
    "    dropoutrate = 0.20\n",
    "    widths = (500, 100, 100, 50)\n",
    "    \n",
    "    # Configure mode\n",
    "    input_dim = 19\n",
    "    \n",
    "    # Architecture definition\n",
    "    model = Sequential()\n",
    "    model.add(Dense(widths[0], input_dim=input_dim, activation=\"selu\"))\n",
    "    model.add(Lambda(lambda x: K.dropout(x, level=dropoutrate)))\n",
    "    for i in range(1,len(widths)):\n",
    "        model.add(Dense(widths[i], activation=\"relu\"))\n",
    "        model.add(Lambda(lambda x: K.dropout(x, level=dropoutrate)))\n",
    "    model.add(Dense(2, activation=\"softmax\"))\n",
    "    \n",
    "    # Compile\n",
    "    model.compile(optimizer=\"adam\", loss=\"binary_crossentropy\", metrics=[\"accuracy\"])\n",
    "    \n",
    "    # Callbacks\n",
    "    callbacks = []\n",
    "    callbacks.append( ModelCheckpoint(model_name+\".hdf5\", save_best_only=True, verbose=verbose) )\n",
    "    callbacks.append( EarlyStopping(patience=10) )\n",
    "    \n",
    "    # Train\n",
    "    history = model.fit(Xtrain, ytrain,\n",
    "                        validation_data=(Xval, yval),\n",
    "                        epochs=epochs,\n",
    "                        callbacks=callbacks,\n",
    "                        verbose=verbose,\n",
    "                        batch_size=1000)\n",
    "\n",
    "    # Load best model\n",
    "    model.load_weights(model_name+\".hdf5\")\n",
    "    \n",
    "    return (model, history)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Test Model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_model(model, Xtest, npasses=100, verbose=0):\n",
    "    results = np.zeros((len(Xtest), npasses))\n",
    "    for i in range(npasses):\n",
    "        results[:,i] = model.predict(Xtest, batch_size=10000)[:,1].flatten()\n",
    "        if verbose > 0:\n",
    "            print(\"Test pass {} succeeded\".format(i))\n",
    "    ypred = np.mean(results, axis=1)\n",
    "    sigma = np.std(1.0*(results>0.5), axis=1)\n",
    "    #sigma = np.std(results, axis=1)\n",
    "    info = np.abs(ypred-0.5)/sigma\n",
    "    #info = np.amax(info) - info\n",
    "    return (ypred, info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Show model uncertainty results**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def uncertainty_results(truth, prediction, info_measure, x=0, y=2):\n",
    "    print(\"{:<20}{}\".format(\"truth.shape\", ytest.shape))\n",
    "    print(\"{:<20}{}\".format(\"pred.shape\", ypred.shape))\n",
    "    print(\"{:<20}{}\".format(\"info.shape\", ysigma.shape))\n",
    "\n",
    "    f, a = plt.subplots(2,2, figsize=(16,16))\n",
    "    a[0,0].scatter(Xtest[:,x], Xtest[:,y], c=ytest[:,1].ravel(), cmap=\"seismic_r\")\n",
    "    a[0,0].set_title(\"Truth\")\n",
    "    a[0,1].scatter(Xtest[:,x], Xtest[:,y], c=ypred.ravel(), cmap=\"seismic_r\")\n",
    "    a[0,1].set_title(\"Prediction\")\n",
    "    a[1,0].scatter(Xtest[:,x], Xtest[:,y], c=np.abs(ytest[:,1]-ypred), cmap=\"Reds\")\n",
    "    a[1,0].set_title(\"Difference\")\n",
    "    a[1,1].scatter(Xtest[:,x], Xtest[:,y], c=ysigma, cmap=\"Purples\")\n",
    "    a[1,1].set_title(\"Uncertainty\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Test model training and testing\n",
    "**Generate Data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xtrain, ytrain = generate_data(10000)\n",
    "Xval, yval = generate_data(50000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Train Model**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model, history = train_model(Xtrain, ytrain, Xval, yval, epochs=20, verbose=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Show resulting model and training history**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show resulting model and history of training\n",
    "model.summary()\n",
    "\n",
    "plt.plot(history.epoch, history.history[\"loss\"], label=\"loss\")\n",
    "plt.plot(history.epoch, history.history[\"val_loss\"], label=\"val_loss\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"epoch\")\n",
    "plt.ylabel(\"MSE\")\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "### Active learning\n",
    "**Configuration**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_start = 10000       # Start size\n",
    "size_iter = 2500         # Number of data points added in each step\n",
    "size_sample = 100000     # Size of set to be checked for uncertainty\n",
    "size_max = 100000        # Maximum size of data set\n",
    "size_val = 1000000       # Size of validation set\n",
    "size_test = 1000000      # Size of test set\n",
    "nmembers = 25            # Committee size\n",
    "niterations = 7          # Number of iterations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Active Sampling function**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def active_sampling(model, Nquery, Nselect, Npasses, random_fraction=0.5):\n",
    "    # Get uncertainty measure\n",
    "    X, _ = generate_data(Nquery)\n",
    "    predictions, info = test_model(model, X, Npasses)   \n",
    "    keysort = np.argsort(info)#[::-1]\n",
    "\n",
    "    # Select actively\n",
    "    select_active = round((1-random_fraction)*Nselect)\n",
    "    selected = X[keysort[:select_active]]\n",
    "    method = np.zeros(Nselect)\n",
    "    method[:select_active] = 1.0\n",
    "\n",
    "    # Add random\n",
    "    select_random = Nselect - select_active\n",
    "    X, _ = generate_data(select_random)\n",
    "    selected = np.vstack((selected, X))\n",
    "\n",
    "    # Label and return\n",
    "    prediction = oracle(selected)\n",
    "    return (selected, prediction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Get model performance**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model_performance(model, Xtest, ytest, niterations):\n",
    "    ypred = np.zeros((len(Xtest),2))\n",
    "    for i in range(niterations):\n",
    "        ypred += model.predict(Xtest, batch_size=10000, verbose=0)\n",
    "    ypred /= niterations\n",
    "    bce = log_loss(ytest, ypred)\n",
    "    acc = accuracy_score(ytest, 1.0*(ypred>0.5))\n",
    "    return {\"bce\":bce, \"acc\":acc}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Create logbooks**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_AL = open(\"log_active_learning.csv\", \"w\")\n",
    "log_AL.write(\"iteration,size,bce,acc\\n\")\n",
    "log_AL.flush()\n",
    "\n",
    "log_RS = open(\"log_random_sampling.csv\", \"w\")\n",
    "log_RS.write(\"iteration,size,bce,acc\\n\")\n",
    "log_RS.flush()\n",
    "\n",
    "def log_result(log, iteration, size, bce, acc):\n",
    "    log.write(\"{},{},{},{}\\n\".format(iteration, size, bce, acc))\n",
    "    log.flush()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Run Active Learning**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for iteration in range(niterations):\n",
    "    print(\"ITERATION {}\".format(iteration))\n",
    "\n",
    "    # Data set creation\n",
    "    Xtrain_AL, ytrain_AL = generate_data(size_start)\n",
    "    Xtrain_RS, ytrain_RS = deepcopy(Xtrain_AL), deepcopy(ytrain_AL)\n",
    "    Xval, yval = generate_data(size_val)\n",
    "    Xtest, ytest = generate_data(size_test)\n",
    "    \n",
    "    while len(Xtrain_AL) < size_max:        \n",
    "        \"\"\" ACTIVE LEARNING \"\"\"\n",
    "        # Train model for active learning\n",
    "        model, _ = train_model(Xtrain_AL, ytrain_AL, Xval, yval, epochs=500, verbose=0)\n",
    "        # Test model performance\n",
    "        performance_AL = model_performance(model, Xtest, ytest, nmembers)\n",
    "        # Store performance\n",
    "        log_result(log_AL, iteration, len(Xtrain_AL), performance_AL[\"bce\"], performance_AL[\"acc\"])\n",
    "        # Active sampling of new points\n",
    "        Xnew, ynew = active_sampling(model, size_sample, size_iter, nmembers, random_fraction=0.5)\n",
    "        # Append new points to active learning\n",
    "        Xtrain_AL = np.vstack((Xtrain_AL, Xnew))\n",
    "        ytrain_AL = np.vstack((ytrain_AL, ynew))\n",
    "        \n",
    "        \"\"\" RANDOM SAMPLING \"\"\"\n",
    "        # Train model for random sampling\n",
    "        model, _ = train_model(Xtrain_RS, ytrain_RS, Xval, yval, epochs=500, verbose=0)\n",
    "        # Test model performance\n",
    "        performance_RS = model_performance(model, Xtest, ytest, nmembers)\n",
    "        # Store performance\n",
    "        log_result(log_RS, iteration, len(Xtrain_RS), performance_RS[\"bce\"], performance_RS[\"acc\"])\n",
    "        # Sample new points\n",
    "        Xnew, ynew = generate_data(size_iter)\n",
    "        # Append new points to random sampling\n",
    "        Xtrain_RS = np.vstack((Xtrain_RS, Xnew))\n",
    "        ytrain_RS = np.vstack((ytrain_RS, ynew))\n",
    "        \n",
    "        \"\"\" LOG AND OUTPUT RESULTS \"\"\"\n",
    "        # Screen\n",
    "        print(\"iteration: {:<5}     size: {:<5}    al-bce: {:<10}    al-acc: {:<10}    rs-bce: {:<10}    rs-acc: {:<10}\".format(\n",
    "            iteration,\n",
    "            len(Xtrain_AL)-size_iter,\n",
    "            round(performance_AL[\"bce\"],6),\n",
    "            round(performance_RS[\"acc\"],6),\n",
    "            round(performance_AL[\"bce\"],6),\n",
    "            round(performance_RS[\"acc\"],6)\n",
    "        ))\n",
    "        # Output data sets to file\n",
    "        al = np.hstack((Xtrain_AL, ytrain_AL))\n",
    "        rs = np.hstack((Xtrain_RS, ytrain_RS))\n",
    "        np.savetxt(\"susyai/active_{}.csv\".format(len(Xtrain_AL)), al, delimiter=\",\")\n",
    "        np.savetxt(\"susyai/random_{}.csv\".format(len(Xtrain_RS)), rs, delimiter=\",\")\n",
    "        \n",
    "        \"\"\" Clear Keras Session \"\"\"\n",
    "        K.clear_session()"
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}