MBL_ETH_Transition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import sklearn\n",
    "# suppress tensorflow compilation warnings\n",
    "import os\n",
    "os.environ['KMP_DUPLICATE_LIB_OK']='True'\n",
    "import tensorflow as tf\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'\n",
    "import numpy as np\n",
    "seed = 0\n",
    "np.random.seed(seed) # fix random seed\n",
    "# tf.set_random_seed(seed)\n",
    "import matplotlib.pyplot as plt\n",
    "%config InlineBackend.figure_format ='retina' #plot high-res img "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data loading and pre-processing\n",
    "We can use either entanglement spectrum or the absolute values of wavefunction as the input "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the data by specifying the data type, interaction strengths and number of samples\n",
    "def load_data(**args):\n",
    "    # This is only useful if you run the en_spectrum.py script to generate your own data set\n",
    "    d_type = args['type']\n",
    "    L = int(args['len'])\n",
    "    j1 = float(args['j1'])\n",
    "    j2 = float(args['j2'])\n",
    "    runs = int(para['run_1']) # number of samples generated\n",
    "    loaded = np.load(\"data_set/j_\"+str(j1)+\"/1_en_spectrum_L=\"+str(L)+\"_j=\"+str(j1)+\".npz\")\n",
    "    # for ent_specta, the size is determined by the cut position\n",
    "    # for the wavefunction, it's by the system size\n",
    "    dim = loaded[d_type].shape\n",
    "    x = np.zeros((runs, dim[0], dim[1]), dtype='float64')\n",
    "    \n",
    "    # load data for the phase 1\n",
    "    for i in range(runs):\n",
    "        filename = \"data_set/j_\"+str(j1)+\"/\"+str(i+1)+\"_en_spectrum_L=\"+str(L)+\"_j=\"+str(j1)+\".npz\"\n",
    "        loaded = np.load(filename)\n",
    "        if d_type == 'wave':\n",
    "            wav = loaded[d_type]\n",
    "            x[i] = np.multiply(wav, np.conj(wav)).real\n",
    "        else:\n",
    "            x[i] = loaded[d_type]\n",
    "            \n",
    "    p1 = np.concatenate((x[:])) \n",
    "    \n",
    "    # load data for the phase 2\n",
    "    runs = int(para['run_2'])    \n",
    "    x = np.zeros((runs, dim[0], dim[1]), dtype='float64')\n",
    "    for i in range(runs):\n",
    "        filename = \"data_set/j_\"+str(j2)+\"/\"+str(i+1)+\"_en_spectrum_L=\"+str(L)+\"_j=\"+str(j2)+\".npz\"\n",
    "        loaded = np.load(filename)\n",
    "        if d_type == 'wave':\n",
    "            wav = loaded[d_type]\n",
    "            x[i] = np.multiply(wav, np.conj(wav)).real\n",
    "        else:\n",
    "            x[i] = loaded[d_type]\n",
    "            \n",
    "    p2 = np.concatenate((x[:])) \n",
    "\n",
    "    # cast the original data to single precision to speed up\n",
    "    p1 = p1.astype('float32')\n",
    "    p2 = p2.astype('float32')\n",
    "    \n",
    "    return p1, p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "para = {'type':'ent','len':'12', 'run_1':'35', 'run_2':'32', 'j1':'0.01', 'j2':'5.0'}\n",
    "\n",
    "# load data according to parameters\n",
    "p1, p2 = load_data(**para)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(90200, 128) (86920, 128)\n"
     ]
    }
   ],
   "source": [
    "para = {'type':'ent','len':'12', 'run_1':'20', 'run_2':'21', 'j1':'0.05', 'j2':'1.5'}\n",
    "# load data according to parameters\n",
    "pp1, pp2 = load_data(**para)\n",
    "p1 = np.concatenate((p1, pp1))\n",
    "p2 = np.concatenate((p2, pp2))\n",
    "print(p1.shape, p2.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "num_classes = 2\n",
    "# processing data \n",
    "def data_pipeline(p1, p2, test_size = 0.2, encode = True):\n",
    "    # create labels for the two phases\n",
    "    l1 = np.ones(p1.shape[0],dtype=np.int8)\n",
    "    l2 = np.zeros(p2.shape[0],dtype=np.int8)\n",
    "    # combine labels and the data\n",
    "    ph1 = np.column_stack((l1, p1))\n",
    "    ph2 = np.column_stack((l2, p2))\n",
    "    # train_test split\n",
    "    dat = np.concatenate((ph1, ph2))\n",
    "    X_train, X_test, Y_train, Y_test = train_test_split(dat[:,1:], dat[:,:1], test_size=test_size, random_state=1)\n",
    "    # convert class vectors to binary class matrices\n",
    "    if encode == True:\n",
    "        Y_train = tf.keras.utils.to_categorical(Y_train, num_classes)\n",
    "        Y_test = tf.keras.utils.to_categorical(Y_test, num_classes)\n",
    "    \n",
    "    return X_train, X_test, Y_train, Y_test\n",
    "\n",
    "# train test split\n",
    "X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.keras import Sequential\n",
    "from tensorflow.keras.layers import Dense, Dropout, Flatten, BatchNormalization\n",
    "\n",
    "def create_model():\n",
    "    # instantiate model\n",
    "    model = Sequential()\n",
    "    # add a normalization layer to improve robustness\n",
    "    model.add(BatchNormalization(axis=-1, input_shape=(X_test.shape[1],)))\n",
    "    # add a dense all-to-all relu layer\n",
    "    model.add(Dense(400, activation='relu'))\n",
    "    # add a dense all-to-all relu layer\n",
    "    model.add(Dense(50, activation='relu'))\n",
    "    # apply dropout with rate 0.4\n",
    "    model.add(Dropout(0.4))\n",
    "    # soft-max layer\n",
    "    model.add(Dense(num_classes, activation='softmax'))\n",
    "    \n",
    "    # compile the model with loss function and optimizer\n",
    "    model.compile(loss=tf.keras.losses.binary_crossentropy,\n",
    "                  optimizer=tf.keras.optimizers.Adam(),\n",
    "                  metrics=['accuracy'])\n",
    "    \n",
    "    return model\n",
    "\n",
    "\n",
    "from tensorflow.keras.callbacks import EarlyStopping\n",
    "#early termination to prevent overfitting\n",
    "early_stopping = EarlyStopping(monitor='val_loss', patience=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# training parameters\n",
    "batch_size = 200\n",
    "epochs = 14\n",
    "\n",
    "# create the deep neural net\n",
    "model_DNN = create_model()\n",
    "\n",
    "# train DNN and store training info in history\n",
    "history = model_DNN.fit(X_train, Y_train,\n",
    "          batch_size=batch_size,\n",
    "          epochs=epochs,\n",
    "          verbose=0,\n",
    "          validation_data=(X_test, Y_test), callbacks=[early_stopping])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate model\n",
    "score = model_DNN.evaluate(X_test, Y_test, verbose=0)\n",
    "\n",
    "# print performance\n",
    "print('Test loss:', score[0])\n",
    "print('Test accuracy:', score[1])\n",
    "# look into training history\n",
    "# summarize history for accuracy\n",
    "plt.plot(history.history['accuracy'])\n",
    "plt.plot(history.history['val_accuracy'])\n",
    "plt.ylabel('model accuracy')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='best')\n",
    "plt.show()\n",
    "\n",
    "# summarize history for loss\n",
    "plt.plot(history.history['loss'])\n",
    "plt.plot(history.history['val_loss'])\n",
    "plt.ylabel('model loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hyperparameters tuning\n",
    "Use [Keras Tuner](https://blog.tensorflow.org/2020/01/hyperparameter-tuning-with-keras-tuner.html) to search for the ideal parameters, including the discrete parameters such number of layers, number neurons in each layer etc., and continuous parameters like the dropout rate and  learning rate which is perphaps the most important one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import kerastuner as kt\n",
    "from kerastuner.tuners import Hyperband, BayesianOptimization\n",
    "\n",
    "\n",
    "def create_model(hp):\n",
    "    # instantiate model\n",
    "    model = Sequential()\n",
    "    # add a normalization layer to improve robustness\n",
    "    model.add(BatchNormalization(axis=-1, input_shape=(X_test.shape[1],)))\n",
    "    # add a dense all-to-all relu layer\n",
    "#     for i in range(hp.Int('num_layers', 2, 6)):\n",
    "#         model.add(layers.Dense(units=hp.Int('units_' + str(i), 32, 512, 32),\n",
    "#                                activation='relu'))\n",
    "    model.add(Dense(units=hp.Int('layer_1_units',\n",
    "                                        min_value=32,\n",
    "                                        max_value=512,\n",
    "                                        step=64),\n",
    "                           activation='relu'))\n",
    "    model.add(Dense(units=hp.Int('layer_2_units',\n",
    "                                        min_value=32,\n",
    "                                        max_value=512,\n",
    "                                        step=64),\n",
    "                           activation='relu'))\n",
    "    model.add(Dropout(\n",
    "      hp.Float('dropout_rate', 0, 0.5, step=0.1, default=0.5)))\n",
    "\n",
    "    # soft-max layer\n",
    "    model.add(Dense(num_classes, activation='softmax'))\n",
    "    \n",
    "    # compile the model with loss function and optimizer    \n",
    "    model.compile(\n",
    "      optimizer=tf.keras.optimizers.Adam(\n",
    "        hp.Float('learning_rate', 1e-4, 1e-1, sampling='log')),\n",
    "      loss='binary_crossentropy',\n",
    "      metrics=['accuracy'])\n",
    "    \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tuner = kt.Hyperband(\n",
    "    create_model,\n",
    "    objective='val_accuracy',\n",
    "    max_epochs=3,\n",
    "    hyperband_iterations=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b_tuner = BayesianOptimization(\n",
    "    create_model,\n",
    "    objective='val_accuracy',\n",
    "    max_trials=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using partial dateset for hp tuning\n",
    "tuner.search(X_train[:20000], Y_train[:20000],\n",
    "             epochs=2,\n",
    "             validation_data=(X_test[:4000], Y_test[:4000]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp = tuner.get_best_hyperparameters(1)[0]\n",
    "print(hp.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Re-train the model with optimized hyper-parameters\n",
    "Now that we've found the (hopefully) best hyper parameters, let's retrain the model to see its performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_model_hp():\n",
    "    # instantiate model\n",
    "    model = Sequential()\n",
    "    # add a normalization layer to improve robustness\n",
    "    model.add(BatchNormalization(axis=-1, input_shape=(X_test.shape[1],)))\n",
    "    # add a dense all-to-all relu layer\n",
    "    model.add(Dense(288, activation='relu'))\n",
    "    # add a dense all-to-all relu layer\n",
    "    model.add(Dense(352, activation='relu'))\n",
    "    model.add(Dropout(0.1))\n",
    "    # soft-max layer\n",
    "    model.add(Dense(num_classes, activation='softmax'))\n",
    "    \n",
    "    # compile the model with loss function and optimizer\n",
    "    model.compile(loss=tf.keras.losses.binary_crossentropy,\n",
    "                  optimizer=tf.keras.optimizers.Adam(learning_rate = 0.003746593),\n",
    "                  metrics=['accuracy'])\n",
    "    \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# training parameters\n",
    "batch_size = 200\n",
    "epochs = 14\n",
    "\n",
    "# create the deep neural net\n",
    "model_hp = create_model_hp()\n",
    "\n",
    "# train DNN and store training info in history\n",
    "history = model_hp.fit(X_train, Y_train,\n",
    "          batch_size=batch_size,\n",
    "          epochs=epochs,\n",
    "          verbose=0,\n",
    "          validation_data=(X_test, Y_test), callbacks=[early_stopping])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate model\n",
    "new_score = model_hp.evaluate(X_test, Y_test, verbose=0)\n",
    "\n",
    "# print performance\n",
    "print('Test loss:', new_score[0])\n",
    "print('Test accuracy:', new_score[1])\n",
    "\n",
    "# look into training history\n",
    "# summarize history for accuracy\n",
    "plt.plot(history.history['accuracy'])\n",
    "plt.plot(history.history['val_accuracy'])\n",
    "plt.ylabel('model accuracy')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='best')\n",
    "plt.show()\n",
    "\n",
    "# summarize history for loss\n",
    "plt.plot(history.history['loss'])\n",
    "plt.plot(history.history['val_loss'])\n",
    "plt.ylabel('model loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see a slight improve of accurary over the old model and a reduction in the gap between the test and training data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load other data set with different $j$ to see how well the network generalizes\n",
    "According to the energy level statistics, the critical point occurs around $j_c\\approx 0.23$.\n",
    "So we use different $j$ values as data set to see how our network works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "para = {'type':'ent','len':'12', 'run_1':'8', 'run_2':'8', 'j1':'0.1', 'j2':'0.5'}\n",
    "\n",
    "# load data according to parameters\n",
    "p1, p2 = load_data(**para)\n",
    "\n",
    "# train test split\n",
    "X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2, 0.99)\n",
    "print(X_train.shape, X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate model\n",
    "score = model_DNN.evaluate(X_test, Y_test, verbose=0)\n",
    "score_hp = model_hp.evaluate(X_test, Y_test, verbose=0)\n",
    "# print performance\n",
    "print('Test loss:', score[0], score_hp[0])\n",
    "print('Test accuracy:', score[1], score_hp[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reverse transform the one-hot matrix to 1-D list\n",
    "_, _, _, Y_test = data_pipeline(p1, p2, 0.99, encode = False)\n",
    "Y_pred = model_hp.predict_classes(X_test, verbose=0)\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "enc = OneHotEncoder(handle_unknown='ignore')\n",
    "enc.fit(Y_test)\n",
    "Y_test = enc.transform(Y_test).toarray()\n",
    "\n",
    "Y_true = enc.inverse_transform(Y_test)\n",
    "Y_true = Y_true.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "cm = confusion_matrix(Y_true, Y_pred)\n",
    "print(cm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to the convention of sci-kit learn, the $c_{0,0}$ component of the matrix represents the true negatives. So the total positive rate is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm[1,1]+cm[1,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "para = {'type':'ent','len':'12', 'run_1':'8', 'run_2':'8', 'j1':'0.15', 'j2':'0.31'}\n",
    "\n",
    "# load data according to parameters\n",
    "p1, p2 = load_data(**para)\n",
    "\n",
    "# train test split\n",
    "X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2, 0.99)\n",
    "print(X_train.shape, X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate model\n",
    "score = model_DNN.evaluate(X_test, Y_test, verbose=0)\n",
    "score_hp = model_hp.evaluate(X_test, Y_test, verbose=0)\n",
    "# print performance\n",
    "print('Test loss:', score[0], score_hp[0])\n",
    "print('Test accuracy:', score[1], score_hp[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using critical data set for identifying critical points\n",
    "According to the energy level statistics, the critical point occurs around $j_c\\approx 0.23$.\n",
    "To benchmark this result, we use $j=0.21$ and $j=0.25$ as data set close to the critical point to see how our network works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "para = {'type':'ent','len':'12', 'run_1':'8', 'run_2':'8', 'j1':'0.21', 'j2':'0.25'}\n",
    "p1, p2 = load_data(**para)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# train test split\n",
    "X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2, 0.99)\n",
    "print(X_train.shape, X_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate model\n",
    "score = model_DNN.evaluate(X_test, Y_test, verbose=0)\n",
    "score_hp = model_hp.evaluate(X_test, Y_test, verbose=0)\n",
    "# print performance\n",
    "print('Test loss:', score[0], score_hp[0])\n",
    "print('Test accuracy:', score[1], score_hp[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}