RMT_Util.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMT Util\n",
    "\n",
    "### TODO:  refactor to allow pyTorch or Keras Models\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.165565Z",
     "start_time": "2018-08-13T04:49:06.705240Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/charlesmartin14/anaconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: compiletime version 3.5 of module 'tensorflow.python.framework.fast_tensor_util' does not match runtime version 3.6\n",
      "  return f(*args, **kwds)\n",
      "/Users/charlesmartin14/.local/lib/python3.6/site-packages/h5py/__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.6.5 |Anaconda custom (64-bit)| (default, Apr 26 2018, 08:42:37) \n",
      "[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]\n",
      "numpy version 1.14.2\n",
      "scipy version 1.0.1\n",
      "tensforflow version 1.5.0-dev20171116\n",
      "keras version 2.1.5\n",
      "sklearn version 0.19.1\n",
      "powerlaw version 1.4.3\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import pickle, time\n",
    "from copy import deepcopy\n",
    "from shutil import copy\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import keras\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense, Flatten\n",
    "from keras.layers.convolutional import Conv2D, MaxPooling2D\n",
    "from keras.layers.normalization import BatchNormalization\n",
    "from keras.callbacks import TensorBoard, EarlyStopping\n",
    "from keras.optimizers import SGD\n",
    "from keras.initializers import Constant\n",
    "import keras.backend as K\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import sklearn\n",
    "from sklearn.decomposition import TruncatedSVD\n",
    "#from sklearn.random_projection import sparse_random_matrix\n",
    "import scipy as sp\n",
    "from scipy.linalg import svd\n",
    "\n",
    "import powerlaw\n",
    "from scipy import optimize\n",
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "\n",
    "print(sys.version)\n",
    "print(\"numpy version {}\".format(np.__version__))\n",
    "print(\"scipy version {}\".format(sp.__version__))\n",
    "print(\"tensforflow version {}\".format(tf.__version__))\n",
    "print(\"keras version {}\".format(keras.__version__))\n",
    "print(\"sklearn version {}\".format(sklearn.__version__))\n",
    "print(\"powerlaw version {}\".format(powerlaw.__version__))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-15T02:19:36.864514Z",
     "start_time": "2018-04-15T02:19:36.861892Z"
    }
   },
   "source": [
    "### MLP3 Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.176265Z",
     "start_time": "2018-08-13T04:49:10.167498Z"
    }
   },
   "outputs": [],
   "source": [
    "def MLP3_model():\n",
    "    model = Sequential()\n",
    "    model.add(Flatten(input_shape=(28, 28, 3)))\n",
    "    model.add(Dense(512, kernel_initializer='glorot_normal',\n",
    "                    bias_initializer=Constant(0.1), activation='relu'))\n",
    "    model.add(Dense(512, kernel_initializer='glorot_normal',\n",
    "                    bias_initializer=Constant(0.1), activation='relu'))\n",
    "    model.add(Dense(512, kernel_initializer='glorot_normal',\n",
    "                    bias_initializer=Constant(0.1), activation='relu'))\n",
    "    model.add(Dense(10, kernel_initializer='glorot_normal',\n",
    "                    bias_initializer=Constant(0.1), activation='softmax'))\n",
    "\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AlexNet Model\n",
    "\n",
    "#### batchnorm = True   only for last 2 layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.207955Z",
     "start_time": "2018-08-13T04:49:10.178087Z"
    }
   },
   "outputs": [],
   "source": [
    "def Alexnet_Model(batchnorm=False):\n",
    "    model = Sequential()\n",
    "    model.add(Conv2D(96, (5, 5), input_shape=(28, 28, 3), kernel_initializer=\n",
    "                     'glorot_normal', bias_initializer=Constant(0.1), padding=\n",
    "                     'same', activation='relu')) \n",
    "    model.add(MaxPooling2D((3, 3), padding='same'))\n",
    "    model.add(BatchNormalization())\n",
    "    model.add(Conv2D(256, (5, 5), kernel_initializer='glorot_normal',\n",
    "                     bias_initializer=Constant(0.1), padding='same',\n",
    "                     activation='relu')) \n",
    "    model.add(MaxPooling2D((3, 3), padding='same'))\n",
    "    model.add(BatchNormalization())\n",
    "    model.add(Flatten())\n",
    "\n",
    "    model.add(Dense(384, kernel_initializer='glorot_normal',\n",
    "                    bias_initializer=Constant(0.1), activation='relu'))\n",
    "    if batchnorm:\n",
    "        model.add(BatchNormalization())\n",
    "    model.add(Dense(192, kernel_initializer='glorot_normal',\n",
    "                    bias_initializer=Constant(0.1), activation='relu'))\n",
    "    if batchnorm:\n",
    "        model.add(BatchNormalization())\n",
    "\n",
    "    model.add(Dense(10, kernel_initializer='glorot_normal',\n",
    "                    bias_initializer=Constant(0.1), activation='softmax'))\n",
    "\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generalized Entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.214065Z",
     "start_time": "2018-08-13T04:49:10.209927Z"
    }
   },
   "outputs": [],
   "source": [
    "def load_keras_matrix(model, weightfile, layer=2):\n",
    "    model.load_weights(weightfile)\n",
    "    return get_keras_matrix(model, layer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.218538Z",
     "start_time": "2018-08-13T04:49:10.215793Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_keras_matrix(model, layer=2):\n",
    "    return model.layers[layer].get_weights()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.223272Z",
     "start_time": "2018-08-13T04:49:10.220460Z"
    }
   },
   "outputs": [],
   "source": [
    "# model is defined above\n",
    "def layer_entropy(model, layer=2):\n",
    "    W = get_keras_matrix(model, layer) \n",
    "    return matrix_entropy(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.235194Z",
     "start_time": "2018-08-13T04:49:10.224935Z"
    }
   },
   "outputs": [],
   "source": [
    "# Trace Normalization\n",
    "def matrix_entropy(W):\n",
    "    W=W/np.trace(W)\n",
    "    m = W.shape[1]\n",
    "    #svd = TruncatedSVD(n_components=m-1, n_iter=7, random_state=10)\n",
    "    #svd.fit(W) \n",
    "    #sv = svd.singular_values_\n",
    "    u, sv, vh = svd(W)\n",
    "    \n",
    "    rank = np.linalg.matrix_rank(W)\n",
    "    p = sv*sv\n",
    "    p = p/np.sum(p)\n",
    "    \n",
    "    if (rank==1):\n",
    "        rank=1.000001\n",
    "    entropy = - np.sum(p*np.log(p)) / np.log(rank) \n",
    "    return entropy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read Training Accuracies\n",
    "\n",
    "We can simply parse these from the output file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.241351Z",
     "start_time": "2018-08-13T04:49:10.236816Z"
    }
   },
   "outputs": [],
   "source": [
    "def read_final_test_accuracy(filename):\n",
    "    cmd = \"tail -1 {} | sed -e 's/\\]//' |  sed -e 's/\\[//' \".format(filename)\n",
    "    !$cmd > tmp.out\n",
    "    result = np.loadtxt(\"tmp.out\",delimiter=\",\")\n",
    "    !rm tmp.out\n",
    "    return result[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.248750Z",
     "start_time": "2018-08-13T04:49:10.243250Z"
    }
   },
   "outputs": [],
   "source": [
    "def read_training_accuracies(filename):\n",
    "    cmd = \"cat {} | grep '50000/' | sed -e 's/.*acc://' \".format(filename)\n",
    "    !$cmd > tmp.out\n",
    "    result = np.loadtxt(\"tmp.out\")\n",
    "    !rm tmp.out\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recompute Test Accuracies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.292149Z",
     "start_time": "2018-08-13T04:49:10.250563Z"
    }
   },
   "outputs": [],
   "source": [
    "def read_test_images(datadir=\"../data\"):\n",
    "\n",
    "    img = tf.placeholder(tf.float32, [28, 28, 3])\n",
    "    norm_image = tf.image.per_image_standardization(img)\n",
    "\n",
    "    print(\"Loading test images...\")\n",
    "    cifar10_test_images = []\n",
    "    cifar10_test_labels = []\n",
    "    datafile = \"{}/cifar-10-batches-py/test_batch\".format(datadir)\n",
    "    \n",
    "    test_file = open(datafile, 'rb')\n",
    "    test_dict = pickle.load(test_file, encoding='bytes')\n",
    "\n",
    "    for image, label in zip(test_dict[b'data'], test_dict[b'labels']):\n",
    "        image_red = np.reshape(image[:1024], (32, 32))[2:-2, 2:-2] / 255.0\n",
    "        image_red = np.reshape(image_red, (28, 28, 1))\n",
    "        image_green = np.reshape(image[1024:2048], (32, 32))[2:-2,\n",
    "                                                                2:-2] / 255.0\n",
    "        image_green = np.reshape(image_green, (28, 28, 1))\n",
    "        image_blue = np.reshape(image[2048:3072], (32, 32))[2:-2, 2:-2] / 255.0\n",
    "        image_blue = np.reshape(image_blue, (28, 28, 1))\n",
    "        image_blue = np.reshape(image_blue, (28, 28, 1))\n",
    "        image = np.concatenate([image_red, image_green, image_blue], axis=-1)\n",
    "        image = norm_image.eval(feed_dict={img:image})\n",
    "        cifar10_test_images.append(image)\n",
    "        label = np.identity(10)[label]\n",
    "        cifar10_test_labels.append(label)\n",
    "    test_file.close()\n",
    "    \n",
    "    return np.array(cifar10_test_images),np.array(cifar10_test_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wigner SemiCircle Plots\n",
    "\n",
    "### <font color='orange'>Not available yet </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Marchenko Pastur FIts and Plots\n",
    "\n",
    "$$\\rho(\\lambda)=\\dfrac{Q}{2\\pi\\sigma_{mp}^{2}}\\dfrac{(\\lambda^{+}-\\lambda)(\\lambda^{-}-\\lambda)}{\\lambda}$$\n",
    "\n",
    "\n",
    "$$\\lambda^{\\pm}=\\sigma_{mp}^{2}\\left(1\\pm\\dfrac{1}{\\sqrt Q}\\right)^{2}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.301248Z",
     "start_time": "2018-08-13T04:49:10.294177Z"
    }
   },
   "outputs": [],
   "source": [
    "def marchenko_pastur_pdf(x_min, x_max, Q, sigma=1):\n",
    "    y=1/Q\n",
    "    x=np.arange(x_min,x_max,0.001)\n",
    "\n",
    "    b=np.power(sigma*(1 + np.sqrt(1/Q)),2) # Largest eigenvalue\n",
    "    a=np.power(sigma*(1 - np.sqrt(1/Q)),2) # Smallest eigenvalue\n",
    "    return x, (1/(2*np.pi*sigma*sigma*x*y))*np.sqrt((b-x)*(x-a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-11T04:51:08.364533Z",
     "start_time": "2018-04-11T04:51:08.350446Z"
    }
   },
   "source": [
    "#### Get sigma from $\\lambda_{max}$ for rectangular matrix\n",
    "\n",
    "\n",
    "$$\\sigma^{2}=\\lambda_{max}\\left(1+\\dfrac{1}{\\sqrt{Q}}\\right)^{-2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.307402Z",
     "start_time": "2018-08-13T04:49:10.303043Z"
    }
   },
   "outputs": [],
   "source": [
    "def calc_sigma(Q, evs):\n",
    "    lmax = np.max(evs)\n",
    "    inv_sqQ = 1.0/np.sqrt(Q)\n",
    "    sigma_2 = lmax/np.square(1+inv_sqQ)\n",
    "    sigma = np.sqrt(sigma_2)\n",
    "    return sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.311541Z",
     "start_time": "2018-08-13T04:49:10.308889Z"
    }
   },
   "outputs": [],
   "source": [
    "def calc_lambda_plus(Q, sigma):\n",
    "    return np.power(sigma*(1 + np.sqrt(1/Q)),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.315898Z",
     "start_time": "2018-08-13T04:49:10.313056Z"
    }
   },
   "outputs": [],
   "source": [
    "def calc_lambda_minus(Q, sigma):\n",
    "    return np.power(sigma*(1 - np.sqrt(1/Q)),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.322410Z",
     "start_time": "2018-08-13T04:49:10.317979Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_Q(W):\n",
    "    if W.shape[1] > W.shape[0]:\n",
    "        M, N = W.shape\n",
    "    else:\n",
    "        N, M = W.shape\n",
    "    Q=N/M\n",
    "    return Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.328426Z",
     "start_time": "2018-08-13T04:49:10.324419Z"
    }
   },
   "outputs": [],
   "source": [
    "def hard_rank(W, sv):\n",
    "    \"\"\"hard rank, using tolerance from numerical recipes, not default scipy tolerance\"\"\"\n",
    "    tol = np.max(sv)*np.finfo(np.float32).eps/(2.0*np.sqrt(np.sum(W.shape)+1))\n",
    "    return np.linalg.matrix_rank(W, tol = tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.335252Z",
     "start_time": "2018-08-13T04:49:10.330148Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_eigenvalues(model, weightfile, layer=7):\n",
    "    \"Read the keras weightfile, get weights for the layer, and compute the eigenvalues(-1)\"\n",
    "    W = load_keras_matrix(model, weightfile, layer)\n",
    "    u, sv, vh = svd(W)\n",
    "    eigenvalues = sv*sv\n",
    "    return eigenvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.354542Z",
     "start_time": "2018-08-13T04:49:10.336785Z"
    }
   },
   "outputs": [],
   "source": [
    "# uss FAST SVD method: notice we miss 1 eigenvalue here...using \n",
    "def get_shuffled_eigenvalues(W, layer=7, num=100):\n",
    "    \"get eigenvalues for this model, but shuffled, num times\"\n",
    "    \n",
    "    print(\"get_shuffled_eigenvalues\")\n",
    "    N, M = W.shape[0], W.shape[1]       \n",
    "\n",
    "    if (N<M):\n",
    "        N, M = W.shape[1], W.shape[0] \n",
    "   \n",
    "    eigenvalues = []\n",
    "    for idx in tqdm(range(num)):\n",
    "        W_shuf = W.flatten()\n",
    "        np.random.shuffle(W_shuf)\n",
    "        W_shuf = W_shuf.reshape([N,M])\n",
    "\n",
    "        u, sv, vh = svd(W_shuf)\n",
    "\n",
    "        eigenvalues.extend(sv*sv)\n",
    "        \n",
    "    evals = (np.array(eigenvalues).flatten())\n",
    "    return evals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit Marchenko Pastur to ESD\n",
    "\n",
    "#### set sigma, or compute sigma from Q "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-28T21:53:11.358461Z",
     "start_time": "2018-08-28T21:53:11.299417Z"
    }
   },
   "outputs": [],
   "source": [
    "def plot_ESD_and_fit(eigenvalues=None, model=None,  weightfile=None, \n",
    "                     layer=2,  Q=1.0, num_spikes=0, sigma=None, \n",
    "                     alpha=0.25, color='blue', skip=False, verbose=True):\n",
    "    \"\"\"Plot histogram of eigenvalues, for Q, and fit Marchenk Pastur.  \n",
    "    If no sigma, calculates from maximum eigenvalue (minus spikes)\n",
    "    Can read keras weights from model if specified.  Does not read PyTorch\"\"\"\n",
    "    if eigenvalues is None:\n",
    "        eigenvalues = get_eigenvalues(model, weightfile, layer)\n",
    "   \n",
    "    plt.hist(eigenvalues, bins=100, alpha=alpha, color=color, density=True, label=r'$\\rho_{emp}(\\lambda)$');\n",
    "\n",
    "    if skip:\n",
    "        return\n",
    "    if not (num_spikes):\n",
    "        num_spikes = 0\n",
    "        \n",
    "    # sort, descending order, minus a few max eigenvalues\n",
    "    evals = np.sort(eigenvalues)[::-1][num_spikes:]\n",
    "\n",
    "    if (sigma is None):\n",
    "        sigma=calc_sigma(Q,evals)\n",
    "        \n",
    "    percent_mass = 100.0*(num_spikes)/len(evals)\n",
    "\n",
    "    ev = np.array(evals)\n",
    "    x_min, x_max =  0, np.max(evals)\n",
    "\n",
    "    x, mp = marchenko_pastur_pdf(x_min, x_max, Q, sigma)\n",
    "    if verbose:\n",
    "        print(\"% spikes outside bulk {0:.2f}\".format(percent_mass))\n",
    "        print(\"% sigma {0:.4f}\".format(sigma))\n",
    "\n",
    "    plt.plot(x,mp, linewidth=1, color = 'r', label=\"MP fit\")\n",
    "    plt.title(\" W{} ESD, MP Sigma={}\".format(layer,sigma));\n",
    "    return sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scree Plots:  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.404671Z",
     "start_time": "2018-08-13T04:49:10.398931Z"
    }
   },
   "outputs": [],
   "source": [
    "# Eigenvalues for correlation weight matrices, for last 2 layers\n",
    "# for Scree plots\n",
    "# check normalization...maybe replace with svd approach above\n",
    "def matrix_eigenvalues(model, layer=2):    \n",
    "    W = model.layers[layer].get_weights()[0]\n",
    "    W = W / np.linalg.norm(W)\n",
    "    WW = np.dot(W.transpose(),W)#/float(N7)\n",
    "    evs, _ = np.linalg.eig(WW)    \n",
    "    return evs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-29T16:24:28.702346Z",
     "start_time": "2018-08-29T16:24:28.682082Z"
    }
   },
   "outputs": [],
   "source": [
    "def scree_plot(model, weightfile, layer=2, color='blue',label=''):    \n",
    "    model.load_weights(weightfile)\n",
    "    evs = matrix_eigenvalues(model, layer)\n",
    "    eigvals = np.flip(np.sort(evs), axis=0)\n",
    "    sing_vals = np.arange(len(eigvals)) + 1\n",
    "    plt.plot(sing_vals, eigvals, color, linewidth=1,label=label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Soft / Stable Rank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.415516Z",
     "start_time": "2018-08-13T04:49:10.412996Z"
    }
   },
   "outputs": [],
   "source": [
    "def stable_rank(evals):\n",
    "    return np.sum(evals)/np.max(evals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.420882Z",
     "start_time": "2018-08-13T04:49:10.417850Z"
    }
   },
   "outputs": [],
   "source": [
    "def layer_soft_rank(model, layer):\n",
    "    W = model.layers[layer].get_weights()[0]\n",
    "    return matrix_soft_rank(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.425551Z",
     "start_time": "2018-08-13T04:49:10.422721Z"
    }
   },
   "outputs": [],
   "source": [
    "def keras_layer_soft_rank(model, layer):\n",
    "    W = keras_matrix(model, layer)\n",
    "    return matrix_soft_rank(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.431082Z",
     "start_time": "2018-08-13T04:49:10.427536Z"
    }
   },
   "outputs": [],
   "source": [
    "def matrix_soft_rank(W):\n",
    "    W=W/np.trace(W)\n",
    "    u, sv, vh = svd(W)\n",
    "    return stable_rank(sv*sv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.439352Z",
     "start_time": "2018-08-13T04:49:10.432733Z"
    }
   },
   "outputs": [],
   "source": [
    "def mp_soft_rank(evals, num_spikes):\n",
    "    evals = np.array(evals)\n",
    "    lambda_max = np.max(evals)\n",
    "    if num_spikes> 0:\n",
    "        evals = np.sort(evals)[::-1][num_spikes:]\n",
    "        lambda_plus = np.max(evals)\n",
    "    else:\n",
    "        lambda_plus = lambda_max\n",
    "        \n",
    "    return lambda_plus/lambda_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.445449Z",
     "start_time": "2018-08-13T04:49:10.441166Z"
    }
   },
   "outputs": [],
   "source": [
    "def calc_mp_soft_rank(evals, Q, sigma):\n",
    "    \n",
    "    lambda_plus = calc_lambda_plus(Q,sigma)\n",
    "    lambda_max = np.max(evals)\n",
    "  \n",
    "    return lambda_plus/lambda_max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-22T20:06:54.614799Z",
     "start_time": "2018-05-22T20:06:54.609018Z"
    }
   },
   "source": [
    "### Eigenvector Localization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.450118Z",
     "start_time": "2018-08-13T04:49:10.446873Z"
    }
   },
   "outputs": [],
   "source": [
    "def keras_matrix(model, modelfile, layer):\n",
    "    model.load_weights(modelfile)\n",
    "    W = model.layers[layer].get_weights()[0]\n",
    "    return W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.454639Z",
     "start_time": "2018-08-13T04:49:10.451722Z"
    }
   },
   "outputs": [],
   "source": [
    "def pytorch_matrix(model, modelfile, layer):\n",
    "   # model.load_weights(modelfile)\n",
    "   # W = model.layers[layer].get_weights()[0]\n",
    "   # return W\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.465013Z",
     "start_time": "2018-08-13T04:49:10.456232Z"
    }
   },
   "outputs": [],
   "source": [
    "def eigenspectrum(W, norm=False):\n",
    "    N, M = W.shape\n",
    "    if N < M:\n",
    "        X = np.dot(W,W.transpose())\n",
    "        \n",
    "    else:\n",
    "        X = np.dot(W.transpose(),W)\n",
    "       \n",
    "    # assumes matrix is full rank\n",
    "    if norm:\n",
    "        X = X/np.trace(X)\n",
    "        \n",
    "    return np.linalg.eig(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.473524Z",
     "start_time": "2018-08-13T04:49:10.466920Z"
    }
   },
   "outputs": [],
   "source": [
    "# TODO: replace with full svd\n",
    "def singular_spectrum(W, norm=False): \n",
    "    if norm:\n",
    "        W = W/np.trace(W)\n",
    "    M = np.min(W.shape)\n",
    "        \n",
    "    svd = TruncatedSVD(n_components=M-1, n_iter=7, random_state=10)\n",
    "    svd.fit(W) \n",
    "    svals = svd.singular_values_\n",
    "    svecs = svd.components_\n",
    "    return svals, svecs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TODO:  fix ensemble runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.483571Z",
     "start_time": "2018-08-13T04:49:10.475281Z"
    }
   },
   "outputs": [],
   "source": [
    "def keras_ensemble_spectrum(model, ensemble_weightfile, layer, num_runs):\n",
    "    \"\"\"compute list of arrays of eigenvals and eigenves for an ensemle run\"\"\"\n",
    "    \n",
    "    eigenvalues, eigenvectors = [], []\n",
    "    for i in tqdm(range(1,num_runs+1)):\n",
    "        filename = ensemble_weightfile.format(i)\n",
    "        W = keras_matrix(model, filename, layer)\n",
    "        evs, evecs = eigenspectrum(W)\n",
    "        \n",
    "        eigenvalues.append(evs)\n",
    "        eigenvectors.append(evecs)\n",
    "    return eigenvalues, eigenvectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.493657Z",
     "start_time": "2018-08-13T04:49:10.485047Z"
    }
   },
   "outputs": [],
   "source": [
    "def keras_singular_spectrum(model, ensemble_weightfile, layer, num_runs):\n",
    "    \"\"\"compute list of arrays of singular vals and vectors for an ensemle run\"\"\"\n",
    "    \n",
    "    svalues, svectors = [], []\n",
    "    for i in tqdm(range(1,num_runs+1)):\n",
    "        filename = ensemble_weightfile.format(i)\n",
    "        W = keras_matrix(model, filename, layer)\n",
    "        svs, svecs = singular_spectrum(W)\n",
    "        \n",
    "        svalues.append(svs)\n",
    "        svectors.append(svecs)\n",
    "    return svalues, svectors "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.498937Z",
     "start_time": "2018-08-13T04:49:10.495980Z"
    }
   },
   "outputs": [],
   "source": [
    "def localization_ratio(vec):\n",
    "    return np.linalg.norm(vec, ord=1) / np.linalg.norm(vec, ord=np.inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.503297Z",
     "start_time": "2018-08-13T04:49:10.500574Z"
    }
   },
   "outputs": [],
   "source": [
    "def participation_ratio(vec):\n",
    "    return np.linalg.norm(vec, ord=2) / np.linalg.norm(vec, ord=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.508812Z",
     "start_time": "2018-08-13T04:49:10.504871Z"
    }
   },
   "outputs": [],
   "source": [
    "def vector_entropy(u):\n",
    "    \"\"\"Vector entropy, as in the  Porter-Thomas distribution\n",
    "    \n",
    "    see:  https://arxiv.org/pdf/cond-mat/9810255.pdf\n",
    "    \"\"\"\n",
    "    u2 = np.linalg.norm(u)/2\n",
    "    entropy = np.sum(np.exp(-u2)/np.sqrt(2.0*np.pi))\n",
    "    return entropy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-25T21:34:59.124609Z",
     "start_time": "2018-05-25T21:34:59.120514Z"
    }
   },
   "source": [
    "### TODO:  Add perEpoch Methods\n",
    "\n",
    "- entropies\n",
    "- ranks\n",
    "- recompute test data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.516240Z",
     "start_time": "2018-08-13T04:49:10.510627Z"
    }
   },
   "outputs": [],
   "source": [
    "def entropy_per_epoch(model, per_epoch_weightfile, num_epochs, layer):\n",
    "    \"\"\"Computer layer entropy per epoch for a keras layer\"\"\"\n",
    "    entropies = []\n",
    "\n",
    "    for epoch in tqdm(range(num_epochs)):\n",
    "        model.load_weights(per_epoch_weightfile.format(epoch))  \n",
    "        entropies.append(layer_entropy(model, layer))\n",
    "\n",
    "    return np.array(entropies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.523473Z",
     "start_time": "2018-08-13T04:49:10.517697Z"
    }
   },
   "outputs": [],
   "source": [
    "def stable_rank_per_epoch(model, per_epoch_weightfile, num_epochs, layer):\n",
    "    \"\"\"Computer layer entropy per epoch for a keras layer\"\"\"\n",
    "    ranks = []\n",
    "\n",
    "    for epoch in tqdm(range(num_epochs)):\n",
    "        model.load_weights(per_epoch_weightfile.format(epoch))  \n",
    "        entropies.append(layer_stable_rank(model, layer))\n",
    "\n",
    "    return np.array(ranks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.534732Z",
     "start_time": "2018-08-13T04:49:10.525049Z"
    }
   },
   "outputs": [],
   "source": [
    "def discrete_entropy(vec, num_bins=100):\n",
    "    vec = vec - np.mean(vec)\n",
    "    h = np.histogram(vec, density=True, bins=num_bins)[0];\n",
    "    p = np.array(h)+0.0000000001\n",
    "    \n",
    "    p = p / np.sqrt(2*np.pi)\n",
    "    p = p / np.sum(p)\n",
    "\n",
    "    #p = p/(2*np.pi)\n",
    "    entropy = -np.sum(p*np.log(p))\n",
    "    entropy = entropy#/(2*np.pi)#/float(num_bins)\n",
    "    return entropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.542449Z",
     "start_time": "2018-08-13T04:49:10.536184Z"
    }
   },
   "outputs": [],
   "source": [
    "def max_discrete_entropy(len_vec,num_bins=100, sample_size = 100000):\n",
    "    \"\"\"compute maximum possible entropy for this vector length and bin sizes\"\"\"\n",
    "    entropies = []\n",
    "    for i in tqdm(range(sample_size)):\n",
    "        test_vec = np.random.normal(0,1,len_vec)\n",
    "        s = discrete_entropy(test_vec, num_bins=num_bins)\n",
    "        entropies.append(s)\n",
    "\n",
    "    return np.max(entropies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Misc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.547615Z",
     "start_time": "2018-08-13T04:49:10.544796Z"
    }
   },
   "outputs": [],
   "source": [
    "def spike_min(Q):\n",
    "    \"\"\"minimum perturbation to W to see a spike\"\"\"\n",
    "    return 1/np.sqrt(np.sqrt(Q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.552487Z",
     "start_time": "2018-08-13T04:49:10.549196Z"
    }
   },
   "outputs": [],
   "source": [
    "def spike_lmax(S,Q):\n",
    "    \"\"\"Maximum spike given a perturbation\"\"\"\n",
    "    S2 = S*S\n",
    "    return ((1.0/Q)+S2)*(1+(1.0/S2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit PowerLaw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-04T00:19:24.191760Z",
     "start_time": "2018-08-04T00:19:24.182535Z"
    }
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.558938Z",
     "start_time": "2018-08-13T04:49:10.553905Z"
    }
   },
   "outputs": [],
   "source": [
    "def mu4alpha(alpha):\n",
    "    if (alpha >= 2.0 and alpha <= 3.0):\n",
    "        mu = 2.0\n",
    "    elif (alpha < 2.0):\n",
    "        mu = 2.0*(alpha - 1.0)\n",
    "    elif (alpha > 3.0):\n",
    "        mu = alpha - 1.0\n",
    "    return mu\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.585319Z",
     "start_time": "2018-08-13T04:49:10.560863Z"
    }
   },
   "outputs": [],
   "source": [
    "def best_dist(fit):\n",
    "    distName = 'power_law'\n",
    "    dist = \"PL\"\n",
    "\n",
    "    R, p = fit.distribution_compare('truncated_power_law', 'power_law', normalized_ratio=True)\n",
    "    if R>0 and p <= 0.05:\n",
    "        distName = 'truncated_power_law'\n",
    "        dist = 'TPL'\n",
    "        \n",
    "    R, p = fit.distribution_compare(distName, 'exponential', normalized_ratio=True)\n",
    "    if R<0 and p <= 0.05:\n",
    "        dist = 'EXP'\n",
    "        return dist\n",
    "\n",
    "    R, p = fit.distribution_compare(distName, 'stretched_exponential', normalized_ratio=True)\n",
    "    if R<0 and p <= 0.05:\n",
    "        dist = 'S_EXP'\n",
    "        return dist\n",
    "        \n",
    "    R, p = fit.distribution_compare(distName, 'lognormal', normalized_ratio=True)\n",
    "    if R<0 and p <= 0.05:\n",
    "        dist = 'LOG_N'\n",
    "        return dist\n",
    "\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.590231Z",
     "start_time": "2018-08-13T04:49:10.586976Z"
    }
   },
   "outputs": [],
   "source": [
    "def fit_powerlaw(evals, verbose=True):\n",
    "    fit = powerlaw.Fit(evals, xmax=np.max(evals))\n",
    "    return [fit.alpha, fit.D, best_dist(fit)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.603504Z",
     "start_time": "2018-08-13T04:49:10.591775Z"
    }
   },
   "outputs": [],
   "source": [
    "# TODO: Try all reshapings\n",
    "def reshape_tensor(W):\n",
    "    \"reshape tensor into a maximum rectangle.  Usually the max dimension (N) is the last, so just preserve it\"\n",
    "    dims = W.shape\n",
    "    N = np.max(dims)\n",
    "    print(dims)\n",
    "\n",
    "    M = 1\n",
    "    for d in dims:\n",
    "        M = M*d\n",
    "    M = int(M / N)\n",
    "    \n",
    "    if dims[-1]==N:\n",
    "        Ws = np.reshape(W,(M,N))\n",
    "    else:\n",
    "        Ws = np.reshape(W,(N,M))\n",
    "        \n",
    "    print(Ws.shape)\n",
    "    return Ws"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_and_plot_powerlaw(evals, verbose=True):\n",
    "    fit = powerlaw.Fit(evals, xmax=np.max(evals))\n",
    "    \n",
    "    alpha, D, best_pl = fit.alpha, fit.D, best_dist(fit)\n",
    "    \n",
    "    print(\"alpha {:3g}, D {:3g}, best_pl  {}\".format(alpha, D, best_pl))\n",
    "    \n",
    "    fig2 = fit.plot_pdf(color='b', linewidth=2)\n",
    "    fit.power_law.plot_pdf(color='b', linestyle='--', ax=fig2)\n",
    "    fit.plot_ccdf(color='r', linewidth=2, ax=fig2)\n",
    "    fit.power_law.plot_ccdf(color='r', linestyle='--', ax=fig2)\n",
    "    plt.show()           \n",
    "            \n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Auto fit MP \n",
    "\n",
    "#### Using Kernel Density Estimator\n",
    "\n",
    "### TODO: float num_spikes, or grid search as int"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.610580Z",
     "start_time": "2018-08-13T04:49:10.605177Z"
    }
   },
   "outputs": [],
   "source": [
    "def marchenko_pastur_fun(x, Q, sigma=1):\n",
    "    y=1/Q\n",
    "    \n",
    "    b=np.power(sigma*(1 + np.sqrt(1/Q)),2) # Largest eigenvalue\n",
    "    a=np.power(sigma*(1 - np.sqrt(1/Q)),2) # Smallest eigenvalue\n",
    "    return x, (1/(2*np.pi*sigma*sigma*x*y))*np.sqrt((b-x)*(x-a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.636830Z",
     "start_time": "2018-08-13T04:49:10.612259Z"
    }
   },
   "outputs": [],
   "source": [
    "def resid_mp(p, evals, Q, num_spikes=0, bw=0.1, debug=False):  \n",
    "    \"residual that floats sigma but NOT Q or num_spikes YET, 10% cutoff each edge\"\n",
    "    sigma = p\n",
    "\n",
    "    # kernel density estimator\n",
    "    kde = KernelDensity(kernel='linear', bandwidth=bw).fit(evals.reshape(-1, 1))\n",
    "    xde =  np.linspace(0, np.max(evals)+0.5, 1000)\n",
    "    X_plot =xde[:, np.newaxis]\n",
    "    log_dens = kde.score_samples(X_plot)\n",
    "    yde = np.exp(log_dens)   \n",
    "    \n",
    "    # MP fit for this sigma\n",
    "    xmp, ymp = marchenko_pastur_fun(xde, Q=Q, sigma=sigma)\n",
    "    \n",
    "    # form residual, remove nan's \n",
    "    resid = ymp-yde\n",
    "    resid = np.nan_to_num(resid)\n",
    "    \n",
    "    if debug:\n",
    "        plt.plot(xde,yde)\n",
    "        plt.plot(xmp,ymp)\n",
    "        plt.show()\n",
    "        print(\"sigma {}  mean residual {}\".format(sigma,np.mean(resid)))\n",
    "\n",
    "    return resid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-13T04:49:10.642577Z",
     "start_time": "2018-08-13T04:49:10.638579Z"
    }
   },
   "outputs": [],
   "source": [
    "def fit_mp(evals, Q):\n",
    "    \"simple fit of evals, only floats sigma right now\"\n",
    "    sigma0=1.0\n",
    "    [sigma1],cov,infodict,mesg,ierr   = optimize.leastsq(resid_mp, [sigma0], args=(evals, Q), full_output=True)\n",
    "    return sigma1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_conv2D_Wmats(Wtensor):\n",
    "    \"\"\"Extract W slices from a 4 index conv2D tensor of shape: (N,M,i,j) or (M,N,i,j).  Return ij (N x M) matrices\"\"\"\n",
    "    Wmats = []\n",
    "    s = Wtensor.shape\n",
    "    N, M, imax, jmax = s[0],s[1],s[2],s[3]\n",
    "    #print(\"tensor shape\", N,M,imax,jmax)\n",
    "    for i in range(imax):\n",
    "        for j in range(jmax):\n",
    "            W = Wtensor[:,:,i,j]\n",
    "            if N < M:\n",
    "                W = W.T\n",
    "            Wmats.append(W)\n",
    "    return Wmats"
   ]
  },
  {
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}