EMF_RBM_EVAL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  EMF RBM Eval\n",
    "\n",
    "evaluate emf_rbm by testing the final class accuracy\n",
    "\n",
    "note: we need the mnist.h5 file, which is the exact same data as can be used in the julia code\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T17:37:29.528287",
     "start_time": "2016-10-20T17:37:27.814344"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/charlesmartin14/anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import h5py\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline \n",
    "\n",
    "from sklearn import linear_model, datasets, metrics, preprocessing \n",
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn.metrics import accuracy_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T17:41:52.628751",
     "start_time": "2016-10-20T17:41:52.622219"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'divide': 'warn', 'invalid': 'warn', 'over': 'warn', 'under': 'warn'}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.datasets import fetch_mldata\n",
    "from sklearn.utils.validation import assert_all_finite\n",
    "from scipy.sparse import csc_matrix, csr_matrix, lil_matrix\n",
    "from sklearn.preprocessing import Binarizer\n",
    "np.seterr(all='warn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T17:58:17.500511",
     "start_time": "2016-10-20T17:58:14.716463"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('List of arrays in this file: \\n', [u'HDF5.name___X', u'HDF5.name___y'])\n",
      "(60000, 784) (60000,)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2117.6342254841179"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hf =  h5py.File('mnist.h5','r')\n",
    "print('List of arrays in this file: \\n', hf.keys())\n",
    "X = np.array(hf.get('HDF5.name___X'))\n",
    "y = np.array(hf.get('HDF5.name___y'))\n",
    "\n",
    "print X.shape, y.shape\n",
    "hf.close()\n",
    "np.linalg.norm(X, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T17:58:17.505566",
     "start_time": "2016-10-20T17:58:17.501895"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.utils.fixes import expit    \n",
    "from sklearn.utils.extmath import safe_sparse_dot\n",
    "\n",
    "def sig_means(x, b, W):\n",
    "    a = safe_sparse_dot(x, W.T) + b\n",
    "    return expit(a, out=a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T17:56:14.156132",
     "start_time": "2016-10-20T17:56:14.153682"
    }
   },
   "source": [
    "### Baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T18:00:54.114142",
     "start_time": "2016-10-20T17:58:21.181275"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BernoulliRBM(batch_size=10, learning_rate=0.1, n_components=256, n_iter=10,\n",
       "       random_state=None, verbose=0)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.neural_network import BernoulliRBM\n",
    "B_rbm = BernoulliRBM()\n",
    "B_rbm.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T18:27:51.192061",
     "start_time": "2016-10-20T18:22:41.668302"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5000 0.945833333333\n"
     ]
    }
   ],
   "source": [
    "p = sig_means(X,B_rbm.intercept_hidden_,B_rbm.components_)\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(p, y, test_size=0.2,random_state=0)\n",
    "for c in [5000]:\n",
    "    lr  = linear_model.LogisticRegression()\n",
    "    lr.C = c\n",
    "    lr.fit(X_train, Y_train)\n",
    "    Y_test_pred = lr.predict(X_test)\n",
    "    acc = accuracy_score(Y_test, Y_test_pred)\n",
    "\n",
    "    print c, acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T18:09:08.723678",
     "start_time": "2016-10-20T18:00:54.115882"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5000 0.945833333333\n"
     ]
    }
   ],
   "source": [
    "p = sig_means(X,B_rbm.intercept_hidden_,B_rbm.components_)\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(p, y, test_size=0.2,random_state=0)\n",
    "for c in [5000]:\n",
    "    lr  = linear_model.LogisticRegression()\n",
    "    lr.C = c\n",
    "    lr.fit(X_train, Y_train)\n",
    "    Y_test_pred = lr.predict(X_test)\n",
    "    acc = accuracy_score(Y_test, Y_test_pred)\n",
    "\n",
    "    print c, acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T18:10:00.754700",
     "start_time": "2016-10-20T18:09:59.421970"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# %load emf_rbm.py\n",
    "\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "\n",
    "\n",
    "from sklearn.base import BaseEstimator\n",
    "from sklearn.base import TransformerMixin\n",
    "from sklearn.externals.six.moves import xrange\n",
    "from sklearn.utils import check_array\n",
    "from sklearn.utils import check_random_state\n",
    "from sklearn.utils import gen_even_slices\n",
    "from sklearn.utils import issparse\n",
    "from sklearn.utils.validation import check_is_fitted\n",
    "\n",
    "from sklearn.utils.fixes import expit  # logistic function  \n",
    "from sklearn.utils.extmath import safe_sparse_dot, log_logistic, softmax\n",
    "\n",
    "class EMF_RBM(BaseEstimator, TransformerMixin):\n",
    "    \"\"\"Extended Mean Field Restricted Boltzmann Machine (RBM).\n",
    "    A Restricted Boltzmann Machine with binary visible units and\n",
    "    binary hidden units. Parameters are estimated using the Extended Mean\n",
    "    Field model, based on the TAP equations\n",
    "    Read more in the :ref:`User Guide <rbm>`.\n",
    "    Parameters\n",
    "    ----------\n",
    "    n_components : int, optional\n",
    "        Number of binary hidden units.\n",
    "    learning_rate : float, optional\n",
    "        The learning rate for weight updates. It is *highly* recommended\n",
    "        to tune this hyper-parameter. Reasonable values are in the\n",
    "        10**[0., -3.] range.\n",
    "    batch_size : int, optional\n",
    "        Number of examples per minibatch.\n",
    "    momentum : float, optional\n",
    "        gradient momentum parameter\n",
    "    decay : float, optional\n",
    "        decay for weight update regularizer\n",
    "    weight_decay: string, optional []'L1', 'L2', None]\n",
    "        weight update regularizer\n",
    "\n",
    "    neq_steps: int, optional\n",
    "        Number of equilibration steps\n",
    "    n_iter : int, optional\n",
    "        Number of iterations/sweeps over the training dataset to perform\n",
    "        during training.\n",
    "    sigma: float, optional\n",
    "        variance of initial W weight matrix\n",
    "    thresh: float, optional\n",
    "        threshold for values in W weight matrix, vectors\n",
    "    verbose : int, optional\n",
    "        The verbosity level. The default, zero, means silent mode.\n",
    "    random_state : integer or numpy.RandomState, optional\n",
    "        A random number generator instance to define the state of the\n",
    "        random permutations generator. If an integer is given, it fixes the\n",
    "        seed. Defaults to the global numpy random number generator.\n",
    "    Attributes\n",
    "    ----------\n",
    "    h_bias : array-like, shape (n_components,)\n",
    "        Biases of the hidden units.\n",
    "    v_bias : array-like, shape (n_features,)\n",
    "        Biases of the visible units.\n",
    "    W : array-like, shape (n_components, n_features)\n",
    "        Weight matrix, where n_features in the number of\n",
    "        visible units and n_components is the number of hidden units.\n",
    "    Examples\n",
    "    --------\n",
    "    >>> import numpy as np\n",
    "    >>> X = np.array([[0, 0, 0], [0, 1, 1], [1, 0, 1], [1, 1, 1]])\n",
    "    >>> model = EMF_RBM(n_components=2)\n",
    "    >>> model.fit(X)\n",
    "    EmfRBM(batch_size=10, learning_rate=0.1, n_components=2, n_iter=10,\n",
    "           random_state=None, verbose=0)\n",
    "    References\n",
    "    ----------\n",
    "    [1] Marylou Gabrie´, Eric W. Tramel1 and Florent Krzakala1, \n",
    "        Training Restricted Boltzmann Machines via the Thouless-Anderson-Palmer Free Energy\n",
    "        https://arxiv.org/pdf/1506.02914\n",
    "    \"\"\"\n",
    "    def __init__(self, n_components=256, learning_rate=0.005, batch_size=100, sigma=0.001, neq_steps = 3,\n",
    "                 n_iter=20, verbose=0, random_state=None, momentum = 0.5, decay = 0.01, weight_decay='L1', thresh=1e-8):\n",
    "        self.n_components = n_components\n",
    "        self.learning_rate = learning_rate\n",
    "        self.batch_size = batch_size\n",
    "        self.n_iter = n_iter\n",
    "        self.verbose = verbose\n",
    "\n",
    "        self.momentum = momentum\n",
    "        self.decay = decay\n",
    "        self.weight_decay = weight_decay\n",
    "\n",
    "        self.sigma = sigma\n",
    "        self.neq_steps = neq_steps\n",
    "\n",
    "        # learning rate / mini_batch\n",
    "        self.lr = learning_rate\n",
    "\n",
    "        # threshold for floats\n",
    "        self.thresh = thresh\n",
    "\n",
    "        # store in case we want to reset\n",
    "        self.random_state = random_state\n",
    "        \n",
    "\n",
    "        # self.random_state_ = random_state\n",
    "        # always start with new random state\n",
    "        self.random_state = check_random_state(random_state)\n",
    "        \n",
    "        # h bias\n",
    "        self.h_bias = np.zeros(self.n_components, )\n",
    "        self.h_samples_ = np.zeros((self.batch_size, self.n_components))\n",
    "        # moved to fit\n",
    "        \n",
    "        self.W = None\n",
    "        self.dW_prev = None\n",
    "        self.W2 = None\n",
    "        self.v_bias = None\n",
    "        \n",
    "\n",
    "    def init_weights(self, X):\n",
    "        \"\"\" If the user specifies the training dataset, it can be useful to                                                                                   \n",
    "        initialize the visibile biases according to the empirical expected                                                                                \n",
    "        feature values of the training data.                                                                                                              \n",
    "\n",
    "        TODO: Generalize this biasing. Currently, the biasing is only written for                                                                         \n",
    "               the case of binary RBMs.\n",
    "        \"\"\"\n",
    "        # \n",
    "        eps = self.thresh\n",
    "\n",
    "        # Mean across  samples \n",
    "        if issparse(X):\n",
    "            probVis = csr_matrix.mean(X, axis=0)\n",
    "        else:\n",
    "            probVis = np.mean(X,axis=0)            \n",
    "\n",
    "        # safe for CSR / sparse mats ?\n",
    "        # do we need it if we use softmax ?\n",
    "        probVis[probVis < eps] = eps            # Some regularization (avoid Inf/NaN)  \n",
    "        #probVis[probVis < (1.0-eps)] = (1.0-eps)   \n",
    "        self.v_bias = np.log(probVis / (1.0-probVis)) # Biasing as the log-proportion\n",
    "        \n",
    "        # (does not work)\n",
    "        # self.v_bias = softmax(probVis)\n",
    "        \n",
    "        # initialize arrays to 0\n",
    "        self.W = np.asarray(\n",
    "            self.random_state.normal(\n",
    "                0,\n",
    "                self.sigma,\n",
    "                (self.n_components, X.shape[1])\n",
    "            ),\n",
    "            order='fortran')\n",
    "\n",
    "        self.dW_prev = np.zeros_like(self.W)\n",
    "        self.W2 = self.W*self.W\n",
    "        return 0\n",
    "\n",
    "\n",
    "    def sample_layer(self, layer):\n",
    "        \"\"\"Sample from the conditional distribution P(h|v) or P(v|h)\"\"\"\n",
    "        self.random_state = check_random_state(self.random_state)\n",
    "        sample = (self.random_state.random_sample(size=layer.shape) < layer) \n",
    "        return sample\n",
    "\n",
    "    def _sample_hiddens(self, v):\n",
    "        \"\"\"Sample from the conditional distribution P(h|v).\n",
    "        Parameters\n",
    "        ----------\n",
    "        v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer to sample from.\n",
    "        Returns\n",
    "        -------\n",
    "        h : array-like, shape (n_samples, n_components)\n",
    "            Values of the hidden layer.\n",
    "        \"\"\"\n",
    "        return self.sample_layer(self._mean_hiddens(v))\n",
    "\n",
    "    def _mean_hiddens(self, v):\n",
    "        \"\"\"Computes the conditional probabilities P(h=1|v).\n",
    "        Parameters\n",
    "        ----------\n",
    "        v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer.\n",
    "        Returns\n",
    "        -------\n",
    "        h : array-like, shape (n_samples, n_components)\n",
    "            Corresponding mean field values for the hidden layer.\n",
    "        \"\"\"\n",
    "        p = safe_sparse_dot(v, self.W.T) + self.h_bias\n",
    "        return expit(p, out=p)\n",
    "\n",
    "    def _sample_visibles(self, h):\n",
    "        \"\"\"Sample from the distribution P(v|h).\n",
    "        Parameters\n",
    "        ----------\n",
    "        h : array-like, shape (n_samples, n_components)\n",
    "            Values of the hidden layer to sample from.\n",
    "        Returns\n",
    "        -------\n",
    "        v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer.\n",
    "        \"\"\"\n",
    "        return sample_layer(self._mean_visible(h))\n",
    "\n",
    "    def _mean_visibles(self, h):\n",
    "        \"\"\"Computes the conditional probabilities P(v=1|h).\n",
    "        Parameters\n",
    "        ----------\n",
    "        h : array-like, shape (n_samples, n_components)\n",
    "            Corresponding mean field values for the hidden layer.\n",
    "        Returns\n",
    "        -------\n",
    "         v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer.     \n",
    "        \"\"\"\n",
    "        #p = np.dot(h, self.W) + self.v_bias\n",
    "        p = safe_sparse_dot(h, W) + self.v_bias\n",
    "        return expit(p, out=p)\n",
    "\n",
    "    def sigma_means(self, x, b, W):\n",
    "        \"\"\"helper class for computing Wx+b \"\"\"\n",
    "        a = safe_sparse_dot(x, W.T) + b\n",
    "        return expit(a, out=a)\n",
    "\n",
    "    def init_batch(self, vis):\n",
    "        \"\"\"initialize the batch for EMF only\"\"\"\n",
    "        v_pos = vis\n",
    "        v_init = v_pos\n",
    "\n",
    "        h_pos = self._mean_hiddens(v_pos)\n",
    "        h_init = h_pos\n",
    "\n",
    "        return v_pos, h_pos, v_init, h_init\n",
    "\n",
    "    def equilibrate(self, v0, h0, iters=3):\n",
    "        \"\"\"Run iters steps of the TAP fixed point equations\"\"\"\n",
    "        mv = v0\n",
    "        mh = h0\n",
    "     \n",
    "        for i in range(iters):\n",
    "            mv = 0.5 *self.mv_update(mv, mh) + 0.5*mv\n",
    "            mh = 0.5 *self.mh_update(mv, mh) + 0.5*mh\n",
    "        return mv, mh\n",
    "\n",
    "    def mv_update(self, v, h):  \n",
    "        \"\"\"update TAP visbile magnetizations, to second order\"\"\"\n",
    "        \n",
    "        # a = np.dot(h, self.W) + self.v_bias\n",
    "        a = safe_sparse_dot(h, self.W) + self.v_bias\n",
    "\n",
    "        h_fluc = h-np.multiply(h,h)\n",
    "        #a += h_fluc.dot(self.W2)*(0.5-v)\n",
    "        \n",
    "        # 0.5-v is elementwise => dense\n",
    "        if issparse(v):\n",
    "            v_half = (0.5-v.todense())\n",
    "        else:\n",
    "             v_half = (0.5-v)\n",
    "            \n",
    "        a += np.multiply(safe_sparse_dot(h_fluc,self.W2), v_half)\n",
    "        return expit(a, out=a)\n",
    "\n",
    "    def mh_update(self, v, h):\n",
    "        \"\"\"update TAP hidden magnetizations, to second order\"\"\"\n",
    "        a = safe_sparse_dot(v, self.W.T) + self.h_bias\n",
    " \n",
    "        v_fluc = (v-(np.multiply(v,v)))\n",
    "        #a += (v-v*v).dot((self.W2).T)*(0.5-h)\n",
    "        \n",
    "        if issparse(h):\n",
    "            h_half = (0.5-h.to_dense())\n",
    "        else:        \n",
    "            h_half = (0.5-h)\n",
    "            \n",
    "        a += np.multiply(safe_sparse_dot(v_fluc,self.W2.T),h_half)\n",
    "\n",
    "        return expit(a, out=a)\n",
    "\n",
    "\n",
    "    def weight_gradient(self, v_pos, h_pos ,v_neg, h_neg):\n",
    "        \"\"\"compute weight gradient of the TAP Free Energy, to second order\"\"\"\n",
    "        # naive  / mean field\n",
    "        dW = safe_sparse_dot(v_pos.T, h_pos, dense_output=True).T - np.dot(h_neg.T, v_neg)\n",
    "        \n",
    "        # tap2 correction\n",
    "        #  elementwise multiplies\n",
    "        h_fluc = (h_neg - np.multiply(h_neg,h_neg)).T\n",
    "        v_fluc = (v_neg - np.multiply(v_neg,v_neg))\n",
    "        #  dW_tap2 = h_fluc.dot(v_fluc)*self.W\n",
    "        dW_tap2 = np.multiply(safe_sparse_dot(h_fluc,v_fluc),self.W)\n",
    "\n",
    "        dW -= dW_tap2\n",
    "        return dW\n",
    "\n",
    "    def score_samples(self, X):\n",
    "        \"\"\"Compute the pseudo-likelihood of X.\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : {array-like, sparse matrix} shape (n_samples, n_features)\n",
    "            Values of the visible layer. Must be all-boolean (not checked).\n",
    "        Returns\n",
    "        -------\n",
    "        pseudo_likelihood : array-like, shape (n_samples,)\n",
    "            Value of the pseudo-likelihood (proxy for likelihood).\n",
    "        Notes\n",
    "        -----\n",
    "        This method is not deterministic: it computes the TAP Free Energy on X,\n",
    "        then on a randomly corrupted version of X, and\n",
    "        returns the log of the logistic function of the difference.\n",
    "        \"\"\"\n",
    "        check_is_fitted(self, \"W\")\n",
    "\n",
    "        v = check_array(X, accept_sparse='csr')\n",
    "        v, v_ = self._corrupt_data(v)       \n",
    "\n",
    "        fe = self._free_energy(v)\n",
    "        fe_ = self._free_energy(v_)\n",
    "        return v.shape[1] * log_logistic(fe_ - fe)\n",
    "    \n",
    "    \n",
    "    def score_samples_TAP(self, X):\n",
    "        \"\"\"Compute the pseudo-likelihood of X using second order TAP\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : {array-like, sparse matrix} shape (n_samples, n_features)\n",
    "            Values of the visible layer. Must be all-boolean (not checked).\n",
    "        Returns\n",
    "        -------\n",
    "        pseudo_likelihood : array-like, shape (n_samples,)\n",
    "            Value of the pseudo-likelihood (proxy for likelihood).\n",
    "        Notes\n",
    "        -----\n",
    "        This method is not deterministic: it computes the TAP Free Energy on X,\n",
    "        then on a randomly corrupted version of X, and\n",
    "        returns the log of the logistic function of the difference.\n",
    "        \"\"\"\n",
    "        check_is_fitted(self, \"W\")\n",
    "\n",
    "        v = check_array(X, accept_sparse='csr')      \n",
    "        v, v_ = self._corrupt_data(v)       \n",
    "\n",
    "        fe = self._free_energy_TAP(v)\n",
    "        fe_ = self._free_energy_TAP(v_)\n",
    "        return v.shape[1] * log_logistic(fe_ - fe)\n",
    "    \n",
    "    def _corrupt_data(self, v):\n",
    "        self.random_state = check_random_state(self.random_state)\n",
    "        \"\"\"Randomly corrupt one feature in each sample in v.\"\"\"\n",
    "        ind = (np.arange(v.shape[0]),\n",
    "               self.random_state.randint(0, v.shape[1], v.shape[0]))\n",
    "        if issparse(v):\n",
    "            data = -2 * v[ind] + 1\n",
    "            v_ = v + sp.csr_matrix((data.A.ravel(), ind), shape=v.shape)\n",
    "        else:\n",
    "            v_ = v.copy()\n",
    "            v_[ind] = 1 - v_[ind]\n",
    "        return v, v_\n",
    "    \n",
    "    \n",
    "    def score_samples_entropy(self, X):\n",
    "        \"\"\"Compute the entropy of X\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : {array-like, sparse matrix} shape (n_samples, n_features)\n",
    "            Values of the visible layer. Must be all-boolean (not checked).\n",
    "        Returns\n",
    "        -------\n",
    "        entropy : array-like, shape (n_samples,)\n",
    "            Value of the entropy.\n",
    "        Notes\n",
    "        -----\n",
    "        This method is not deterministic: it computes the entropy on X,\n",
    "        then on a randomly corrupted version of X, and returns the difference.\n",
    "        \"\"\"\n",
    "        check_is_fitted(self, \"W\")\n",
    "\n",
    "        v = check_array(X, accept_sparse='csr')\n",
    "        v, v_ = self._corrupt_data(v)       \n",
    "\n",
    "        s = self._entropy(v)\n",
    "        s_ = self._entropy(v_)\n",
    "        return v.shape[1] * (s_ - s)\n",
    "\n",
    "    \n",
    "        #TODO: run per column\n",
    "    def _denoise(self, m, eps=1.0e-8):\n",
    "        \"\"\"denoise magnetization\"\"\"\n",
    "        m = np.maximum(m,eps)\n",
    "        m = np.minimum(m,1.0-eps)\n",
    "        return m\n",
    "\n",
    "\n",
    "    def _free_energy_TAP(self, v):\n",
    "        \"\"\"Computes the TAP Free Energy F(v) to second order\n",
    "        Parameters\n",
    "        ----------\n",
    "        v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer.\n",
    "        Returns\n",
    "        -------\n",
    "        free_energy : array-like, shape (n_samples,)\n",
    "            The value of the free energy.\n",
    "        \"\"\"\n",
    "        fe = (- safe_sparse_dot(v, self.v_bias)\n",
    "                - np.logaddexp(0, safe_sparse_dot(v, self.W.T)\n",
    "                               + self.h_bias).sum(axis=1))\n",
    "        \n",
    "        h = self._mean_hiddens(v)\n",
    "        mv, mh = self.equilibrate(v, h, iters=self.neq_steps)\n",
    "        \n",
    "        mv = self._denoise(mv)\n",
    "        mh = self._denoise(mh)\n",
    "\n",
    "        # sum over nodes: axis=1\n",
    "        \n",
    "        U_naive = (-safe_sparse_dot(mv, self.v_bias) \n",
    "                    -safe_sparse_dot(mh, self.h_bias) \n",
    "                        -(mv.dot(self.W.T)*(mh)).sum(axis=1))     \n",
    "\n",
    "        Entropy = ( -(mv*np.log(mv)+(1.0-mv)*np.log(1.0-mv)).sum(axis=1)  \n",
    "                    -(mh*np.log(mh)+(1.0-mh)*np.log(1.0-mh)).sum(axis=1) )\n",
    "                   \n",
    "        h_fluc = (mh - (mh*mh))\n",
    "        v_fluc = (mv - (mv*mv))\n",
    "\n",
    "        # if we do it this way, we need to normalize by 1/batch_size \n",
    "        # which we need to obtain from the W2 matrix\n",
    "        # (I think because of the double sum)\n",
    "        # this is not obvious in the paper...have to be very careful here\n",
    "        tap_norm = 1.0/float(mv.shape[0])\n",
    "        dW_tap2 = h_fluc.dot(self.W2).dot(v_fluc.T)\n",
    "\n",
    "        # julia way, does not require extra norm, but maybe slower ?\n",
    "        # dW_tap2 = h_fluc.dot(self.W2)*v_fluc\n",
    "        Onsager = -0.5*(dW_tap2).sum(axis=1)*tap_norm\n",
    "        fe_tap = U_naive + Onsager - Entropy\n",
    "\n",
    "        return fe_tap - fe\n",
    "\n",
    "\n",
    "    \n",
    "    def _free_energy(self, v):\n",
    "        \"\"\"Computes the RBM Free Energy F(v) \n",
    "        Parameters\n",
    "        ----------\n",
    "        v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer.\n",
    "        Returns\n",
    "        -------\n",
    "        free_energy : array-like, shape (n_samples,)\n",
    "            The value of the free energy.\n",
    "        \"\"\"\n",
    "        fe = (- safe_sparse_dot(v, self.v_bias)\n",
    "                - np.logaddexp(0, safe_sparse_dot(v, self.W.T)\n",
    "                               + self.h_bias).sum(axis=1) )\n",
    "\n",
    "        return fe \n",
    "    \n",
    "    \n",
    "    def _entropy(self, v):\n",
    "        \"\"\"Computes the TAP Free Energy F(v) to second order\n",
    "        Parameters\n",
    "        ----------\n",
    "        v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer.\n",
    "        Returns\n",
    "        -------\n",
    "        entropy : array-like, shape (n_samples,)\n",
    "            The value of the entropy.\n",
    "        \"\"\"\n",
    "         \n",
    "        h = self._mean_hiddens(v)\n",
    "        mv, mh = self.equilibrate(v, h, iters=self.neq_steps)\n",
    "\n",
    "        mv = self._denoise(mv)\n",
    "        mh = self._denoise(mh)\n",
    "\n",
    "        # appears to be wrong ?  unsure why ?  maybe because it is not denoised !!!\n",
    "        Entropy = ( -(mv*np.log(mv)+(1.0-mv)*np.log(1.0-mv)).sum(axis=1)  \n",
    "                    -(mh*np.log(mh)+(1.0-mh)*np.log(1.0-mh)).sum(axis=1)  )\n",
    "                         \n",
    "        return Entropy\n",
    "\n",
    "\n",
    "    \n",
    "    def _free_energy(self, v):\n",
    "        \"\"\"Computes the RBM Free Energy F(v) \n",
    "        Parameters\n",
    "        ----------\n",
    "        v : array-like, shape (n_samples, n_features)\n",
    "            Values of the visible layer.\n",
    "        Returns\n",
    "        -------\n",
    "        free_energy : array-like, shape (n_samples,)\n",
    "            The value of the free energy.\n",
    "        \"\"\"\n",
    "        fe = (- safe_sparse_dot(v, self.v_bias)\n",
    "                - np.logaddexp(0, safe_sparse_dot(v, self.W.T)\n",
    "                               + self.h_bias).sum(axis=1) )\n",
    "\n",
    "        return fe \n",
    "\n",
    "    \n",
    "    def partial_fit(self, X, y=None):\n",
    "        \"\"\"Fit the model to the data X which should contain a partial\n",
    "        segment of the data.\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : array-like, shape (n_samples, n_features)\n",
    "            Training data.\n",
    "        Returns\n",
    "        -------\n",
    "        self : EMF_RBM\n",
    "            The fitted model.\n",
    "        \"\"\"\n",
    "        X = check_array(X, accept_sparse='csr', dtype=np.float64)\n",
    "        if not hasattr(self, 'random_state_'):\n",
    "            self.random_state_ = check_random_state(self.random_state)\n",
    "        if not hasattr(self, 'W'):\n",
    "            self.W = np.asarray(\n",
    "                self.random_state_.normal(\n",
    "                    0,\n",
    "                    0.01,\n",
    "                    (self.n_components, X.shape[1])\n",
    "                ),\n",
    "                order='F')\n",
    "        if not hasattr(self, 'h_bias'):\n",
    "            self.h_bias = np.zeros(self.n_components, )\n",
    "        if not hasattr(self, 'v_bias'):\n",
    "            self.v_bias = np.zeros(X.shape[1], )\n",
    "\n",
    "        # not used ?\n",
    "        #if not hasattr(self, 'h_samples_'):\n",
    "        #    self.h_samples_ = np.zeros((self.batch_size, self.n_components))\n",
    "\n",
    "        self._fit(X)\n",
    "\n",
    "    def _fit(self, v_pos):\n",
    "        \"\"\"Inner fit for one mini-batch.\n",
    "        Adjust the parameters to maximize the likelihood of v using\n",
    "        Extended Mean Field theory (second order TAP equations).\n",
    "        Parameters\n",
    "        ----------\n",
    "        v_pos : array-like, shape (n_samples, n_features)\n",
    "            The data to use for training.\n",
    "        \"\"\"\n",
    "        X_batch = v_pos\n",
    "        lr = float(self.learning_rate) / X_batch.shape[0]\n",
    "        decay = self.decay\n",
    "\n",
    "        v_pos, h_pos, v_init, h_init = self.init_batch(X_batch)\n",
    "              \n",
    "        a = safe_sparse_dot(h_init, self.W, dense_output=True) + self.v_bias\n",
    "        a = expit(a, out=a)\n",
    "\n",
    "        # get_negative_samples\n",
    "        v_neg, h_neg = self.equilibrate(v_init, h_init, iters=self.neq_steps) \n",
    "        \n",
    "        # basic gradient\n",
    "        dW = self.weight_gradient(v_pos, h_pos ,v_neg, h_neg) \n",
    "\n",
    "        # regularization based on weight decay\n",
    "        #  similar to momentum >\n",
    "        if self.weight_decay == \"L1\":\n",
    "            dW -= decay * np.sign(self.W)\n",
    "        elif self.weight_decay == \"L2\":\n",
    "            dW -= decay * self.W\n",
    "\n",
    "        # can we use BLAS here ?\n",
    "        # momentum\n",
    "        # note:  what do we do if lr changes per step ? not ready yet\n",
    "        dW += self.momentum * self.dW_prev  \n",
    "        # update\n",
    "        self.W += lr * dW \n",
    "\n",
    "        # storage for next iteration\n",
    "\n",
    "        # is this is a memory killer \n",
    "        self.dW_prev =  dW  \n",
    "        \n",
    "        # is this wasteful...can we avoid storing 2X the W mat ?\n",
    "        # elementwise multiply\n",
    "        self.W2 = np.multiply(self.W,self.W)\n",
    "\n",
    "        # update bias terms\n",
    "        #   csr matrix sum is screwy, returns [[1,self.n_components]] 2-d array  \n",
    "        #   so I always use np.asarray(X.sum(axis=0)).squeeze()\n",
    "        #   although (I think) this could be optimized\n",
    "        self.v_bias += lr * (np.asarray(v_pos.sum(axis=0)).squeeze() - np.asarray(v_neg.sum(axis=0)).squeeze())\n",
    "        self.h_bias += lr * (np.asarray(h_pos.sum(axis=0)).squeeze() - np.asarray(h_neg.sum(axis=0)).squeeze())\n",
    "\n",
    "        return 0\n",
    "\n",
    "    \n",
    "    def fit(self, X, y=None):\n",
    "        \"\"\"Fit the model to the data X.\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : {array-like, sparse matrix} shape (n_samples, n_features)\n",
    "            Training data.\n",
    "        Returns\n",
    "        -------\n",
    "        self : EMF_RBM\n",
    "            The fitted model.\n",
    "        \"\"\"\n",
    "        verbose = self.verbose\n",
    "        X = check_array(X, accept_sparse='csr', dtype=np.float64)\n",
    "        self.random_state = check_random_state(self.random_state)\n",
    "        \n",
    "        self.init_weights(X)\n",
    "        \n",
    "        n_samples = X.shape[0]\n",
    "        n_batches = int(np.ceil(float(n_samples) / self.batch_size))\n",
    "        \n",
    "\n",
    "        batch_slices = list(gen_even_slices(n_batches * self.batch_size,\n",
    "                                            n_batches, n_samples))\n",
    "        \n",
    "        begin = time.time()\n",
    "        for iteration in xrange(1, self.n_iter + 1):\n",
    "            #print \"iter \", iteration\n",
    "            for batch_slice in batch_slices:\n",
    "                self._fit(X[batch_slice])\n",
    "\n",
    "            #print \"batches done\"\n",
    "            if verbose:\n",
    "                end = time.time()\n",
    "                print(\"[%s] Iteration %d, pseudo-likelihood = %.2f,\"\n",
    "                      \" time = %.2fs\"\n",
    "                      % (type(self).__name__, iteration,\n",
    "                         self.score_samples(X).mean(), end - begin))\n",
    "                begin = end\n",
    "\n",
    "        return self\n",
    "    \n",
    "    def transform(self, X):\n",
    "        \"\"\"Compute the hidden layer activation probabilities, P(h=1|v=X).\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : {array-like, sparse matrix} shape (n_samples, n_features)\n",
    "            The data to be transformed.\n",
    "        Returns\n",
    "        -------\n",
    "        h : array, shape (n_samples, n_components)\n",
    "            Latent representations of the data.\n",
    "        \"\"\"\n",
    "        check_is_fitted(self, \"W\")\n",
    "\n",
    "        X = check_array(X, accept_sparse='csr', dtype=np.float64)\n",
    "        return self._mean_hiddens(X)\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T18:19:30.989710",
     "start_time": "2016-10-20T18:14:45.432509"
    },
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[EMF_RBM] Iteration 1, pseudo-likelihood = -186.39, time = 10.48s\n",
      "[EMF_RBM] Iteration 2, pseudo-likelihood = -150.08, time = 13.96s\n",
      "[EMF_RBM] Iteration 3, pseudo-likelihood = -134.27, time = 13.64s\n",
      "[EMF_RBM] Iteration 4, pseudo-likelihood = -122.06, time = 14.72s\n",
      "[EMF_RBM] Iteration 5, pseudo-likelihood = -115.19, time = 14.35s\n",
      "[EMF_RBM] Iteration 6, pseudo-likelihood = -107.51, time = 14.28s\n",
      "[EMF_RBM] Iteration 7, pseudo-likelihood = -101.65, time = 14.27s\n",
      "[EMF_RBM] Iteration 8, pseudo-likelihood = -98.46, time = 14.36s\n",
      "[EMF_RBM] Iteration 9, pseudo-likelihood = -94.03, time = 14.37s\n",
      "[EMF_RBM] Iteration 10, pseudo-likelihood = -94.02, time = 14.79s\n",
      "[EMF_RBM] Iteration 11, pseudo-likelihood = -90.95, time = 13.80s\n",
      "[EMF_RBM] Iteration 12, pseudo-likelihood = -88.69, time = 14.17s\n",
      "[EMF_RBM] Iteration 13, pseudo-likelihood = -84.68, time = 16.27s\n",
      "[EMF_RBM] Iteration 14, pseudo-likelihood = -84.79, time = 14.33s\n",
      "[EMF_RBM] Iteration 15, pseudo-likelihood = -85.27, time = 14.01s\n",
      "[EMF_RBM] Iteration 16, pseudo-likelihood = -81.95, time = 14.18s\n",
      "[EMF_RBM] Iteration 17, pseudo-likelihood = -80.62, time = 14.97s\n",
      "[EMF_RBM] Iteration 18, pseudo-likelihood = -81.72, time = 14.47s\n",
      "[EMF_RBM] Iteration 19, pseudo-likelihood = -77.75, time = 14.03s\n",
      "[EMF_RBM] Iteration 20, pseudo-likelihood = -78.29, time = 13.98s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "EMF_RBM(batch_size=100, decay=0.01, learning_rate=0.005, momentum=0.5,\n",
       "    n_components=256, n_iter=20, neq_steps=3,\n",
       "    random_state=<mtrand.RandomState object at 0x10ba715a0>, sigma=0.001,\n",
       "    thresh=1e-08, verbose=True, weight_decay='L1')"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "emf_rbm = EMF_RBM(verbose=True)\n",
    "emf_rbm.fit(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-20T18:27:51.192061",
     "start_time": "2016-10-20T18:22:41.668302"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5000 0.95075\n"
     ]
    }
   ],
   "source": [
    "p = sig_means(X,emf_rbm.h_bias,emf_rbm.W)\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(p, y, test_size=0.2,random_state=0)\n",
    "for c in [5000]:\n",
    "    lr  = linear_model.LogisticRegression()\n",
    "    lr.C = c\n",
    "    lr.fit(X_train, Y_train)\n",
    "    Y_test_pred = lr.predict(X_test)\n",
    "    acc = accuracy_score(Y_test, Y_test_pred)\n",
    "\n",
    "    print c, acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}