EMF_RBM_Test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EMF RBM Class Test\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:42:55.357201",
     "start_time": "2016-10-30T18:42:55.350473"
    },
    "collapsed": false
   },
   "outputs": [],
   "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": 121,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T19:07:11.831734",
     "start_time": "2016-10-30T19:07:10.487863"
    },
    "collapsed": false,
    "scrolled": false
   },
   "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",
    "        print \"mv, mh\", np.linalg.norm(mv,ord=2), np.linalg.norm(mh,ord=2)\n",
    "\n",
    "        \n",
    "        mv = self._denoise(mv)\n",
    "        mh = self._denoise(mh)\n",
    "        \n",
    "        print \"denoised mv, mh\", np.linalg.norm(mv,ord=2), np.linalg.norm(mh,ord=2)\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",
    "        print \"S \", np.mean(Entropy)\n",
    "        print \"U \", np.mean(U_naive)\n",
    "        print \"O \", np.mean(Onsager)\n",
    "\n",
    "        return fe_tap \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",
    "            for batch_slice in batch_slices:\n",
    "                self._fit(X[batch_slice])\n",
    "\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": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-03T20:35:48.241019",
     "start_time": "2016-10-03T20:35:48.236692"
    }
   },
   "source": [
    "# EMF Tests\n",
    "\n",
    "Designed to match results of julia code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-03T20:32:59.404629",
     "start_time": "2016-10-03T20:32:59.402070"
    }
   },
   "source": [
    "## Test EMF class init\n",
    "\n",
    "Xdigits\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:41.957209",
     "start_time": "2016-10-30T18:25:41.820690"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1797, 64)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.datasets import load_digits\n",
    "from sklearn.utils.validation import assert_all_finite\n",
    "from scipy.sparse import csc_matrix, csr_matrix, lil_matrix\n",
    "from sklearn.utils.testing import (assert_almost_equal, assert_array_equal,\n",
    "                                   assert_true)\n",
    "\n",
    "\n",
    "from sklearn.preprocessing import Binarizer\n",
    "np.seterr(all='warn')\n",
    "\n",
    "Xdigits = load_digits().data\n",
    "Xdigits -= Xdigits.min()\n",
    "Xdigits /= Xdigits.max()\n",
    "\n",
    "b = Binarizer(threshold=0.001, copy=True)\n",
    "Xdigits = b.fit_transform(Xdigits)\n",
    "print Xdigits.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### create dataset for julia"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-03T21:08:43.534724",
     "start_time": "2016-10-03T21:08:43.529474"
    }
   },
   "source": [
    "X = Xdigits.copy()\n",
    "hf =  h5py.File('xdigits.h5',\"w\")\n",
    "hf.create_dataset(\"X\", data=X) \n",
    "hf.close()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-03T21:08:50.521176",
     "start_time": "2016-10-03T21:08:50.513828"
    }
   },
   "source": [
    "hf =  h5py.File('xdigits.h5',\"r\")\n",
    "print(\"keys\",hf.keys())\n",
    "X = np.array(hf.get('X'))\n",
    "hf.close()\n",
    "print \"norm of X \",np.linalg.norm(X,ord=2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### julia  xdigits_ex.jl   results checked"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-03T21:37:07.818848",
     "start_time": "2016-10-03T21:37:07.813802"
    }
   },
   "source": [
    "### test init"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:43.818148",
     "start_time": "2016-10-30T18:25:43.795520"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def test_init():\n",
    "    X = Xdigits.copy()\n",
    "    assert_almost_equal(np.linalg.norm(X,ord=2), 211.4983270228649  , decimal=12)\n",
    "\n",
    "    rbm = EMF_RBM(momentum=0.5, n_components=64, batch_size=50 , decay=0.01, learning_rate=0.005, n_iter=0, sigma=0.001, neq_steps=3, verbose=True)\n",
    "    rbm.fit(X)\n",
    "    assert_true(np.linalg.norm(rbm.h_bias, ord=2)==0.0)\n",
    "    assert_true(np.linalg.norm(rbm.lr)==0.005)\n",
    "    assert_true(np.linalg.norm(rbm.momentum)==0.5)\n",
    "    assert_true(np.linalg.norm(rbm.decay)==0.01)\n",
    "    assert_true(np.linalg.norm(rbm.n_iter)==0)\n",
    "    assert_true(np.linalg.norm(rbm.neq_steps)==3)\n",
    "    assert_true(np.linalg.norm(rbm.sigma)==0.001)\n",
    "    assert_true(np.linalg.norm(rbm.verbose)==True)\n",
    "    assert_true(np.linalg.norm(rbm.n_components)==64)\n",
    "    assert_true(np.linalg.norm(rbm.thresh)==1e-8)\n",
    "    assert_true(np.linalg.norm(rbm.batch_size)==50)\n",
    "\n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2), 38.97455, decimal=5)\n",
    "\n",
    "    #assert_true(np.linalg.norm(rbm.weight_decay)=='L1')\n",
    "    assert_array_equal(X, Xdigits)\n",
    "    \n",
    "test_init()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test one batch exactly \n",
    "\n",
    "#### $\\sigma\\sim0$, decay = 0 (no regularization)\n",
    "\n",
    "- note: julia code must also use 1e-8 as eps, not 1e-6\n",
    "\n",
    "norm of W  0.007628825568441182  \n",
    "norm W2: 1.8196753262900838e-6\n",
    "\n",
    "hb 1.004859173557616e-18  \n",
    "vb 38.97452180357592\n",
    "\n",
    "\n",
    "pseudo l-hood: -12.962946404073032  \n",
    "entropy: 68.52758766764042  \n",
    "TAP free_energy: -24.38378040255592  \n",
    "U naive: -48.5436547091488  \n",
    "free energy: -9268.746331749979  \n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:44.545754",
     "start_time": "2016-10-30T18:25:44.485958"
    },
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "\nArrays are not almost equal to 7 decimals\n ACTUAL: -117.07124372005572\n DESIRED: -24.383780402555928",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-67-9ccf3d08276c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     31\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mrbm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     32\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 33\u001b[0;31m \u001b[0mrbm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtest_one_batch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-67-9ccf3d08276c>\u001b[0m in \u001b[0;36mtest_one_batch\u001b[0;34m()\u001b[0m\n\u001b[1;32m     27\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     28\u001b[0m     \u001b[0massert_almost_equal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mavg_entropy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m68.52758766764042\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 29\u001b[0;31m     \u001b[0massert_almost_equal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mavg_free_energy_tap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m24.383780402555928\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     30\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mrbm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/charlesmartin14/anaconda/lib/python2.7/site-packages/numpy/testing/utils.pyc\u001b[0m in \u001b[0;36massert_almost_equal\u001b[0;34m(actual, desired, decimal, err_msg, verbose)\u001b[0m\n\u001b[1;32m    535\u001b[0m         \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    536\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mround\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdesired\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mactual\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdecimal\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 537\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mAssertionError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_build_err_msg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    538\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    539\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: \nArrays are not almost equal to 7 decimals\n ACTUAL: -117.07124372005572\n DESIRED: -24.383780402555928"
     ]
    }
   ],
   "source": [
    "def test_one_batch():\n",
    "    X = Xdigits.copy()\n",
    "    rbm = EMF_RBM(momentum=0.5, n_components=64, batch_size=100,\n",
    "                  decay=0.00, learning_rate=0.005, n_iter=0,\n",
    "                  sigma=1e-16, neq_steps=3, verbose=True, weight_decay=None)\n",
    "    rbm.init_weights(X);\n",
    "      \n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2), 0.0)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W2,ord=2), 0.0 )\n",
    "    assert_almost_equal(np.linalg.norm(rbm.dW_prev,ord=2),0.0 )\n",
    "    \n",
    "    X_batch = X[0:100,:]    \n",
    "    assert_almost_equal(np.linalg.norm(X_batch,ord=2), 49.31032989212045)\n",
    "\n",
    "    rbm.partial_fit(X_batch);\n",
    "    \n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2), 0.007628825568441182 )\n",
    "    \n",
    "    scored_free_energy = np.average(rbm.score_samples(X_batch))\n",
    "    avg_free_energy_tap = np.average(rbm._free_energy_TAP(X_batch))\n",
    "    avg_entropy = np.average(np.average(rbm._entropy(X_batch)))\n",
    "        \n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2), 38.9745218036)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.h_bias,ord=2),0.0 )    \n",
    "    assert_almost_equal(np.linalg.norm(rbm.dW_prev,ord=2), 152.57651136882203)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W2,ord=2), 0.000001819675326290)\n",
    "\n",
    "    assert_almost_equal(avg_entropy, 68.52758766764042)\n",
    "    assert_almost_equal(avg_free_energy_tap, -24.383780402555928)\n",
    "\n",
    "    return rbm\n",
    "\n",
    "rbm = test_one_batch()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-08T11:07:28.040581",
     "start_time": "2016-10-08T11:07:28.032305"
    },
    "collapsed": false
   },
   "source": [
    "### Test 2 batches exactly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:45.249709",
     "start_time": "2016-10-30T18:25:45.221251"
    },
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def test_two_batches():\n",
    "    X = Xdigits.copy()\n",
    "    rbm = EMF_RBM(momentum=0.5, n_components=64, batch_size=100,\n",
    "                  decay=0.00, learning_rate=0.005, n_iter=0,\n",
    "                  sigma=1e-16, neq_steps=3, verbose=True, weight_decay=None)\n",
    "    rbm.init_weights(X);\n",
    "    X_batch = X[0:100,:]     \n",
    "    rbm.partial_fit(X_batch);\n",
    "    X_batch = X[100:200,:]  \n",
    "    rbm.partial_fit(X_batch);\n",
    "\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2), 0.015478158879359825 )\n",
    "    \n",
    "    scored_free_energy = np.average(rbm.score_samples(X_batch))\n",
    "    avg_free_energy_tap = np.average(rbm._free_energy_TAP(X_batch))\n",
    "    avg_entropy = np.average(np.average(rbm._entropy(X_batch)))\n",
    "        \n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2), 38.974504602139554)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.h_bias,ord=2),1.0779652694386856e-6)    \n",
    "    assert_almost_equal(np.linalg.norm(rbm.dW_prev,ord=2), 178.06423558738115)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W2,ord=2), 8.120675004806954e-6)\n",
    "    #assert_almost_equal(avg_entropy, )\n",
    "    #assert_almost_equal(avg_free_energy_tap, )\n",
    "\n",
    "    return rbm\n",
    "\n",
    "rbm = test_two_batches()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test one epoch, $\\sigma=0.001$\n",
    "\n",
    "compare:\n",
    "\n",
    "5 julia runs  \n",
    "batch norm of W, hb, vb  0.015177951725370209 6.125160958113443e-5 38.974531344645186  \n",
    "batch norm of W, hb, vb  0.016005072745766846 6.132506125735679e-5 38.974534343561935  \n",
    "batch norm of W, hb, vb  0.015518275427920199 6.143705375221393e-5 38.97453267232916\n",
    "batch norm of W, hb, vb  0.016618832753491925 6.14604623830071e-5 38.97453303623846\n",
    "batch norm of W, hb, vb  0.015643733669880935 6.131198883353152e-5 38.97453109464897\n",
    "\n",
    "10 BernoulliRBM runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:46.445655",
     "start_time": "2016-10-30T18:25:46.324472"
    },
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def test_one_epoch():\n",
    "    X = Xdigits.copy()\n",
    "    rbm = EMF_RBM(momentum=0.5, n_components=64, batch_size=100,\n",
    "                  decay=0.01, learning_rate=0.005, n_iter=1, \n",
    "                  sigma=0.001, neq_steps=3, verbose=False)\n",
    "    rbm.fit(X);\n",
    "    \n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2),38.974531, decimal=4)\n",
    "    # really between 0.015 and 0.0165: hard to test properly with a single statement\n",
    "\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2),0.0165, decimal=2)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.h_bias,ord=2),0.000061, decimal=2)\n",
    "    \n",
    "    # non tap FE totally wrong\n",
    "    # FE ~ -2x.x\n",
    "\n",
    "    scored_free_energy = np.average(rbm.score_samples(X))\n",
    "    \n",
    "    avg_free_energy_tap = np.average(rbm._free_energy_TAP(X))\n",
    "    avg_entropy = np.average(np.average(rbm._entropy(X)))\n",
    "    \n",
    "    #assert_almost_equal(scored_free_energy, -24, decimal=0)  \n",
    "    #assert_almost_equal(avg_free_energy_tap, -25, decimal=0)  \n",
    "    assert_almost_equal(avg_entropy, 68.8, decimal=0)\n",
    "    return rbm\n",
    "\n",
    "for i in range(1):\n",
    "    test_one_epoch()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-03T22:07:51.818339",
     "start_time": "2016-10-03T22:07:51.814845"
    }
   },
   "source": [
    "### test partial fit  (1 iter, 1 batch)\n",
    "\n",
    "<font color='red'>It might be better to remove the regularization to test this</font>\n",
    "\n",
    "running this 2x works ... unsure why ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:54.020264",
     "start_time": "2016-10-30T18:25:53.542741"
    },
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def test_partial_fit():\n",
    "    X = Xdigits.copy()\n",
    "    rbm = EMF_RBM(momentum=0.5, n_components=64, batch_size=100,\n",
    "                  decay=0.01, learning_rate=0.005, n_iter=0, \n",
    "                  sigma=0.000000001, neq_steps=3, verbose=True)\n",
    "    rbm.init_weights(X);\n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2), 38.9745518)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2), 0.000000001)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W2,ord=2), 0.000000001)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.dW_prev,ord=2), 0.000000001)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.h_bias,ord=2), 0.000000001)\n",
    "\n",
    "    X_batch = Xdigits.copy()[0:100]\n",
    "    assert_almost_equal(np.linalg.norm(X_batch,ord=2), 49.3103298921)\n",
    "    rbm.partial_fit(X_batch)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2),0.007629, decimal=4)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2),38.974521, decimal=4)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.h_bias,ord=2),0.0, decimal=3)    \n",
    "    \n",
    "    #there are large variations in dw_prev\n",
    "    assert_almost_equal(np.linalg.norm(rbm.dW_prev,ord=2),152.6, decimal=1)\n",
    "\n",
    "# test stochastically (sometimes will fail due to roundoff error in dw_prev)\n",
    "for i in range(100):\n",
    "    test_partial_fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-04T16:04:44.744505",
     "start_time": "2016-10-04T16:04:44.742373"
    }
   },
   "source": [
    "### Test 2 iterations of the partial fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:55.465556",
     "start_time": "2016-10-30T18:25:54.884728"
    },
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def test_partial_fit_2iters():\n",
    "    X = Xdigits.copy()\n",
    "    rbm = EMF_RBM(momentum=0.5, n_components=64, batch_size=100,\n",
    "                  decay=0.01, learning_rate=0.005, n_iter=0, \n",
    "                  sigma=0.00000001, neq_steps=3, verbose=True)\n",
    "    rbm.init_weights(X);\n",
    "    X_batch = Xdigits.copy()[0:100]\n",
    "    assert_almost_equal(np.linalg.norm(X_batch,ord=2), 49.3103298921)\n",
    "    rbm.partial_fit(X_batch)\n",
    "    \n",
    "    X_batch = Xdigits.copy()[100:200]\n",
    "    assert_almost_equal(np.linalg.norm(X_batch,ord=2), 48.96867960939811)\n",
    "    rbm.partial_fit(X_batch)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2),38.974504602)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.h_bias,ord=2),0.000001, decimal=6)\n",
    "    \n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2),0.0154, decimal=3)\n",
    "    # not correct ?\n",
    "    assert_almost_equal(np.linalg.norm(rbm.dW_prev,ord=2),177.75, decimal=1)\n",
    "   \n",
    "\n",
    "# test stochastically\n",
    "for i in range(100):\n",
    "    test_partial_fit_2iters()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:56.056337",
     "start_time": "2016-10-30T18:25:55.509063"
    },
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def test_fit_xdigits():\n",
    "    X = Xdigits.copy()\n",
    "    rbm = EMF_RBM(momentum=0.5, n_components=64, batch_size=100,\n",
    "                  decay=0.01, learning_rate=0.005, n_iter=20, \n",
    "                  sigma=0.001, neq_steps=3, verbose=False)\n",
    "    rbm.fit(X);\n",
    "    \n",
    "    \n",
    "    \n",
    "    assert_almost_equal(np.linalg.norm(rbm.W,ord=2),0.02, decimal=1)\n",
    "    assert_almost_equal(np.linalg.norm(rbm.v_bias,ord=2),38.9747, decimal=3)\n",
    "    # why is h so different ?\n",
    "    assert_almost_equal(np.linalg.norm(rbm.h_bias,ord=2),0.0012, decimal=2)\n",
    "    return rbm\n",
    "    \n",
    "rbm = test_fit_xdigits()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### test sample hiddens\n",
    "\n",
    "IDK why there is a divide-by-zero error\n",
    "\n",
    "Perhaps the exitp() needs some regularization ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:57.140226",
     "start_time": "2016-10-30T18:25:57.110442"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/charlesmartin14/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:143: RuntimeWarning: divide by zero encountered in divide\n"
     ]
    }
   ],
   "source": [
    "def test_sample_hiddens():\n",
    "    rng = np.random.RandomState(0)\n",
    "    X = Xdigits[:100]\n",
    "    rbm1 = EMF_RBM(n_components=2, batch_size=5,\n",
    "                        n_iter=5, random_state=42)\n",
    "    rbm1.fit(X)\n",
    "\n",
    "    h = rbm1._mean_hiddens(X[0])\n",
    "    hs = np.mean([rbm1._sample_hiddens(X[0]) for i in range(100)], 0)\n",
    "\n",
    "    assert_almost_equal(h, hs, decimal=1)\n",
    "\n",
    "test_sample_hiddens()  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <font color='red'>why divide by zero error ? <font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test Verbose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:58.805702",
     "start_time": "2016-10-30T18:25:58.587613"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.externals.six.moves import cStringIO as StringIO\n",
    "def test_rbm_verbose():\n",
    "    rbm = EMF_RBM(n_iter=2, verbose=10)\n",
    "    old_stdout = sys.stdout\n",
    "    sys.stdout = StringIO()\n",
    "    try:\n",
    "        rbm.fit(Xdigits)\n",
    "    finally:\n",
    "        sys.stdout = old_stdout\n",
    "        \n",
    "test_rbm_verbose()   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-04T01:29:02.736552",
     "start_time": "2016-10-04T01:29:02.734464"
    }
   },
   "source": [
    "### Test Transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:25:59.475610",
     "start_time": "2016-10-30T18:25:59.442538"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def test_transform():\n",
    "    X = Xdigits[:110] # using 100 causes divide by zero error in mean_hiddens()!\n",
    "    rbm1 = EMF_RBM(n_components=16, batch_size=5,\n",
    "                        n_iter=5, random_state=42)\n",
    "    rbm1.fit(X)\n",
    "\n",
    "    Xt1 = rbm1.transform(X)\n",
    "    Xt2 = rbm1._mean_hiddens(X)\n",
    "\n",
    "    assert_array_equal(Xt1, Xt2)\n",
    "    \n",
    "test_transform()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test means_hidden\n",
    "\n",
    "should compare to older RBM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:26:01.231844",
     "start_time": "2016-10-30T18:26:01.220013"
    },
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.02649814]\n",
      " [ 0.02009084]]\n",
      "[-0.00342014 -0.00348868]\n",
      "[ 0.05]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.neural_network import BernoulliRBM\n",
    "rng = np.random.RandomState(42)\n",
    "X = np.array([[0.], [1.]])\n",
    "rbm = BernoulliRBM(n_components=2, batch_size=2,\n",
    "                    n_iter=42, random_state=rng)\n",
    "# you need that much iters\n",
    "rbm.fit(X)\n",
    "#assert_almost_equal(rbm1.W, np.array([[0.02649814], [0.02009084]]), decimal=4)\n",
    "#assert_almost_equal(rbm1.gibbs(X), X)\n",
    "print rbm.components_\n",
    "print rbm.intercept_hidden_\n",
    "print rbm.intercept_visible_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-05T23:42:12.136222",
     "start_time": "2016-10-05T23:42:12.133731"
    }
   },
   "source": [
    "### <font color='red'>Need to work this out by hand </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:26:02.635752",
     "start_time": "2016-10-30T18:26:02.612794"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.5         0.5       ]\n",
      " [ 0.50012314  0.49996445]]\n",
      "[[ 0.00049258]\n",
      " [-0.0001422 ]]\n",
      "[ -7.96246297e-06]\n"
     ]
    }
   ],
   "source": [
    "def test_mean_hiddens():\n",
    "    # Im not entirely sure why this happens, but the hidden units all go to 1/2\n",
    "    # and the h array is (2,2)\n",
    "    # need to do by hand\n",
    "    rng = np.random.RandomState(42)\n",
    "    X = np.array([[0.], [1.]])\n",
    "    rbm = EMF_RBM(n_components=2, batch_size=2,\n",
    "                        n_iter=42, random_state=rng, \n",
    "                        decay = 0.0, weight_decay=None, momentum=0)\n",
    "    rbm.fit(X)\n",
    "    h = rbm._mean_hiddens(X)\n",
    "    assert_true(h.shape==(2,2))\n",
    "    assert_almost_equal(np.linalg.norm(h,ord=2), 1.0, decimal=4)\n",
    "    assert_almost_equal(h[0,0], 0.5, decimal=3)\n",
    "    assert_almost_equal(h[0,1], 0.5, decimal=3)\n",
    "    assert_almost_equal(h[1,0], 0.5, decimal=3)\n",
    "    assert_almost_equal(h[1,1], 0.5, decimal=3)\n",
    "\n",
    "    print h\n",
    "    print rbm.W\n",
    "    print rbm.v_bias\n",
    "    \n",
    "test_mean_hiddens()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test equlibrate\n",
    "\n",
    "???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:26:03.898991",
     "start_time": "2016-10-30T18:26:03.883674"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  5.15573045e-06]\n",
      " [ -1.97855285e-05]]\n"
     ]
    }
   ],
   "source": [
    "def test_fit_equilibrate():\n",
    "    # Equlibrate on the RBM hidden layer should be able to recreate [[0], [1]]\n",
    "    # from the same input\n",
    "    rng = np.random.RandomState(42)\n",
    "    X = np.array([[0.], [1.]])\n",
    "    rbm1 = EMF_RBM(n_components=2, batch_size=2,\n",
    "                        n_iter=42, random_state=rng)\n",
    "    # you need that much iters\n",
    "    rbm1.fit(X)\n",
    "    #assert_almost_equal(rbm1.W, np.array([[0.02649814], [0.02009084]]), decimal=4)\n",
    "    #assert_almost_equal(rbm1.gibbs(X), X)\n",
    "    return rbm1, X\n",
    "\n",
    "rbm, X = test_fit_equilibrate()\n",
    "print rbm.W"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test using sparse CSR matrix\n",
    "\n",
    "### <font color='red'>need to check sparse matrix results </font>\n",
    "\n",
    "Must confirm that sparse matrix multiplies are not mis-interpeted as matrix dot products"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:26:06.354961",
     "start_time": "2016-10-30T18:26:06.299591"
    },
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/charlesmartin14/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:143: RuntimeWarning: divide by zero encountered in divide\n"
     ]
    }
   ],
   "source": [
    "def test_small_sparse():\n",
    "    # EMF_RBM should work on small sparse matrices.\n",
    "    X = csr_matrix(Xdigits[:4])\n",
    "    EMF_RBM().fit(X)       # no exception\n",
    "\n",
    "test_small_sparse()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test Free Energy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:26:07.771153",
     "start_time": "2016-10-30T18:26:07.764182"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot\n",
    "import matplotlib as mpl\n",
    "%matplotlib inline  \n",
    "\n",
    "def show_image(image): \n",
    "    fig = pyplot.figure()\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    imgplot = ax.imshow(image, cmap=mpl.cm.Greys)\n",
    "    imgplot.set_interpolation('nearest')\n",
    "    ax.xaxis.set_ticks_position('top')\n",
    "    ax.yaxis.set_ticks_position('left')\n",
    "    pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Test Free Energy and Entropy Calculations\n",
    "\n",
    "actual julia output from xdigits_FE.jl\n",
    "\n",
    "(info statements added into monitor; I should find a way to add to test code )\n",
    "\n",
    "S [68.92548157440935,68.67917042062827,68.7382937888165,68.6467445638933,68.94092201361534]  \n",
    "FE <1-5>[-90.41392263605844 -98.57232874119751 -96.67160538171822 -99.72457836849503 -89.84668949506056]  \n",
    "FE TAP <1-5>[-117.38187025746029 -117.39051762052955 -117.39024519247155 -117.39408456128287   -117.37959261213285]\n",
    "\n",
    "\n",
    "?? m vis, hid 11.117791807149395  8.936583879865992  \n",
    "?? denoised m vis, hid 11.11779180715007  8.936583879865992  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T19:29:43.666350",
     "start_time": "2016-10-30T19:29:43.512668"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "entropy  [ 68.92500812  68.67868513  68.73780754  68.64628658  68.94040984]\n",
      "free energies, old  [-90.41429712 -98.57310481 -96.67241348 -99.72444489 -89.84836316]\n",
      "mv, mh 210.707428877 169.423346805\n",
      "denoised mv, mh 210.707428877 169.423346805\n",
      "S  68.8462458135\n",
      "U  -48.5352952066\n",
      "O  -2.95285368833e-06\n",
      "TAP free energies  [-117.38220661 -117.39085455 -117.39058211 -117.3944207  -117.37993056]\n"
     ]
    }
   ],
   "source": [
    "def test_free_energy():\n",
    "    X = Xdigits.copy()\n",
    "    rbm = EMF_RBM(n_iter=1, n_components=64,  decay=0.001, sigma=0.0000000000000001, neq_steps=5)\n",
    "    rbm.fit(X) \n",
    "    \n",
    "    s = rbm._entropy(X)\n",
    "    print \"entropy \", s[0:5]\n",
    "        \n",
    "    fe = rbm._free_energy(X)\n",
    "    print \"free energies, old \" , fe[0:5]\n",
    "    \n",
    "    fe_tap = rbm._free_energy_TAP(X)\n",
    "    print \"TAP free energies \", fe_tap[0:5]\n",
    "    \n",
    "    assert_almost_equal(s[0], 68.92548, decimal=3)\n",
    "    assert_almost_equal(s[1], 68.679170, decimal=3) # a bit more off\n",
    "    \n",
    "    assert_almost_equal(fe[0], -90.4139, decimal=3)\n",
    "    assert_almost_equal(fe[1], -98.5723, decimal=2) # a bit more off\n",
    "  \n",
    "    assert_almost_equal(fe_tap[0], -117.3819, decimal=2)\n",
    "    assert_almost_equal(fe_tap[1], -117.3905, decimal=2)    \n",
    "    \n",
    "    return rbm\n",
    "\n",
    "rbm = test_free_energy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T18:28:07.105745",
     "start_time": "2016-10-30T18:28:07.081681"
    }
   },
   "source": [
    "### IDK why the mv norms are different...weird ... will dig into a bit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T19:34:38.154506",
     "start_time": "2016-10-30T19:34:38.009643"
    },
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAAD7CAYAAAC2TgIoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAC0ZJREFUeJzt3V2MXVUZxvHnGVqaltrGiAFC5UvTGC4UMKExxTiGD/mI\ncGWAaEh6gTeaEmMIhsRMvfES5MIbAlRUQEIjigYJTTAxKB+FtlJoK0oDFKENJKSCvQH7enE2tTnM\nMHuYtfbs/c7/l0zmzJmd864zs59Z6+xZZy1HhADkM7HQDQBQB+EGkiLcQFKEG0iKcANJEW4gqQUL\nt+1Lbe+1/aLtmyrXutP2QdvP1axzTL01th+z/YLtXbY3Vq63zPZTtnc09aZq1mtqTtjebvuh2rWa\nei/b/lvzHJ+uXGu17Qds72l+h+sq1lrbPKftzedDxc6XiOj8Q6M/Kv+UdLqkpZJ2Svp8xXoXSDpH\n0nMdPb+TJZ3T3F4p6e81n19TZ0Xz+ThJT0o6v3K970v6laSHOvqZ7pP0yY5q/VzShub2EkmrOqo7\nIel1SZ8p8XgL1XOfL+kfEfFKRLwn6deSrqpVLCIel/R2rcefpt6BiNjZ3H5X0h5Jp1auebi5uUyj\nE7La7CTbayRdLumOWjWmK6sORpq2V0n6SkRslqSIeD8i/l27buMiSS9FxP4SD7ZQ4T5V0rFP4DVV\nPvkXiu0zNBo1PFW5zoTtHZIOSNoaEdsqlrtV0o2q+AdkGiFpq+1ttq+vWOdMSW/Z3twMlW+3vbxi\nvWNdLem+Ug/GBbWKbK+UtEXSDU0PXk1EHImIcyWtkbTO9tk16ti+QtLBZmTi5qML6yPiPI1GDN+1\nfUGlOksknSfpZ029w5J+WKnWUbaXSrpS0gOlHnOhwv0vSacd8/Wa5r40bC/RKNi/jIjfdVW3GUL+\nSdKllUqsl3Sl7X0a9TJfs/2LSrWOiog3ms9vSnpQo5d2NbwmaX9EPNN8vUWjsNd2maRnm+dXxEKF\ne5ukz9k+3fbxkq6RVPuqa5e9jCTdJWl3RNxWu5DtE22vbm4vl3SxpL01akXEzRFxWkScpdHv7bGI\nuK5GrQ/YXtGMgmT7BEmXSHq+Rq2IOChpv+21zV0XStpdo9aYa1VwSC6NhiCdi4j/2v6epEc1+gNz\nZ0TsqVXP9r2SJiV9yvarkqY+uGBSqd56Sd+StKt5HRySbo6IRyqVPEXS3bYnNPp53h8RD1eqtRBO\nkvSg7dDonL0nIh6tWG+jpHuaofI+SRsq1pLtFRpdTPtO0cdtLsEDSIYLakBShBtIinADSRFuICnC\nDSRV7F9hzb8pACyAiPjQHI6iPffHeefK1NRU9XfblKiFYev7uTmfejNhWA4kRbiBpBY83JOTkylr\nYfi6Pl9K1ys2/dR2ZH5tanf5nhOUlv3crH5BDUB/EG4gKcINJNUq3F0uQwygjFkvqDULALyo0YoU\nr2u0iso1EbF37DguqKG3sp+bH/eCWqfLEAMoo024F80yxEAmXFADkmrzrrDWyxBv2rTp6O3JyUlm\nhAELqM0FteM02uvqQklvSHpa0rXjq5VyQQ19lv3cnO6C2qw9d3S8DDGAMphb3hI997BlPzeZWw4s\nIoQbSIpwA0kRbiApwg0kRbiBpAg3kBThBpIi3EBSxbYTAvqs6xmGfZgRR88NJEW4gaQIN5AU4QaS\nItxAUoQbSIpwA0kRbiApwg0kNWu4bd9p+6Dt57poEIAy2vTcmyV9vXZDAJQ1a7gj4nFJb3fQFgAF\n8ZobSKrou8LYTgjoj1abEtg+XdLvI+ILH3EMmxIAjS6zMN9NCdx8ABiINv8Ku1fSXyWttf2q7Q31\nmwVgvtgrrCWG5ZiLIQ3LAQwM4QaSItxAUoQbSIpwA0kRbiApwg0kRbiBpAg3kNRg9wpjxlhZXc8u\n5PdXHz03kBThBpIi3EBShBtIinADSRFuICnCDSRFuIGkCDeQVJsFEtfYfsz2C7Z32d7YRcMAzM+s\nCyTaPlnSyRGx0/ZKSc9Kuioi9o4d1+kCiUxfLIvpp2UNYoHEiDgQETub2+9K2iPp1PJNBFDSnF5z\n2z5D0jmSnqrRGADltH5XWDMk3yLphqYH/xD2CgP6o+1eYUsk/UHSHyPithmO4TX3gPGau6xBvOZu\n3CVp90zBBtA/ba6Wr5f0Z0m7JEXzcXNEPDJ2HD33gNFzl9WHnnuwe4VlPzm6RrjL6kO4maEGJEW4\ngaQIN5AU4QaSItxAUoQbSIpwA0kRbiApwg0kNdi9wrLresYY8qHnBpIi3EBShBtIinADSRFuICnC\nDSRFuIGkCDeQFOEGkpp1hprtZRotkHh8c/yWiPhx7YYBmJ+265aviIjDto+T9BdJGyPi6bFjWCCx\noOzTT/n9lTOvBRIj4nBzc5lGvXfuMw9IoFW4bU/Y3iHpgKStEbGtbrMAzFerd4VFxBFJ59peJem3\nts+OiN3jx7FXGNAfc96UwPaPJP0nIm4Zu5/X3AXxmnvYBvGa2/aJtlc3t5dLuljS3vJNBFBSm2H5\nKZLutj2h0R+D+yPi4brNAjBf7BXWUwzLh20Qw3IAw0S4gaQIN5AU4QaSItxAUoQbSIpwA0kRbiAp\nwg0kRbiBpAg3kBThBpIi3EBShBtIinADSRFuICnCDSRFuIGkCDeQVOtwNxsTbLf9UM0GAShjLj33\nDZI+tBEBgH5qu53QGkmXS7qjbnMAlNK2575V0o1iA0BgMNrsz32FpIMRsdP2pKQZF5xmrzCgP2bd\nlMD2TyR9W9L7kpZL+oSk30TEdWPHsSlBQWxKMGx92JRgTjuO2P6qpB9ExJXTfI9wF0S4h60P4eb/\n3EBS7BXWU/Tcw0bPDaAawg0kRbiBpAg3kBThBpIi3EBShBtIinADSRFuIKlZ3xWGxSH7jLHFiJ4b\nSIpwA0kRbiApwg0kRbiBpAg3kBThBpIi3EBShBtIqtUMNdsvSzok6Yik9yLi/JqNAjB/baefHpE0\nGRFv12wMgHLaDss9h2MB9EDbwIakrba32b6+ZoMAlNF2WL4+It6w/WmNQr4nIh4fP4i9woD+mPOm\nBLanJL0TEbeM3c+mBAV1vSlB9p9n1waxKYHtFbZXNrdPkHSJpOfLNxFASW2G5SdJetB2NMffExGP\n1m0WgPlir7CeYlg+bIMYlgMYJsINJEW4gaQIN5AU4QaSItxAUoQbSIpwA0kRbiAp9grrKWaMldX1\njL8+oOcGkiLcQFKEG0iKcANJEW4gKcINJEW4gaQIN5AU4QaSahVu26ttP2B7j+0XbK+r3TAA89N2\n+ultkh6OiG/aXiJpRcU2AShg1tVPba+StCMiPjvLcax+it7KPLd8PqufninpLdubbW+3fbvt5eWb\nCKCkNj33lyQ9KenLEfGM7Z9KOhQRU2PHxdTU/++qvVcYPTfmYjH23G3CfZKkJyLirObrCyTdFBHf\nGDuOYTl6azGGe9ZheUQclLTf9trmrgsl7S7cPgCFtdpOyPYXJd0haamkfZI2RMShsWPoudFbi7Hn\nZq8wLAqLMdzMUAOSItxAUoQbSIpwA0kRbiApwg0kRbiBpAg3kBThBpIa7F5hXc84YkZcWZlnjPUF\nPTeQFOEGkiLcQFKEG0iKcANJEW4gKcINJEW4gaRmDbfttbZ3NGuW77B9yPbGLhoH4OOb0xpqtick\nvSZpXUTsH/tep2uodY0ZamVlPle6VmoNtYskvTQebAD9M9dwXy3pvhoNAVBW62G57aWSXpd0dkS8\nOc33GZajtcznStdmGpbP5V1hl0l6drpgf2DTpk1Hb9feKwzAR5tLz32fpEci4u4Zvk/PjdYynytd\nm9eOI7ZXSHpF0lkR8c4MxxButJb5XOlauu2Euka4y8p8rnSN7YSARYZwA0kRbiApwg0kRbiBpAg3\nkBThBpIi3EBShBtIinADSRXdK4wpmmiLc6W+YnPLAfQLw3IgKcINJEW4gaQIN5AU4QaS+h9PbpxX\nSfOeGAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11305f790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "v = X[0,:]\n",
    "show_image(np.reshape(v,[8,8]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TODO:  fix APIs so they can run on 1 vector at a time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T19:34:54.719723",
     "start_time": "2016-10-30T19:34:54.697803"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mv, mh 4.9744512468 3.99674446312\n",
      "denoised mv, mh 4.9744512468 3.99674446312\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "'axis' entry is out of bounds",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-139-6f6eb4191f6d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mrbm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_free_energy_TAP\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-121-fb4089d53c4b>\u001b[0m in \u001b[0;36m_free_energy_TAP\u001b[0;34m(self, v)\u001b[0m\n\u001b[1;32m    424\u001b[0m         U_naive = (-safe_sparse_dot(mv, self.v_bias) \n\u001b[1;32m    425\u001b[0m                     \u001b[0;34m-\u001b[0m\u001b[0msafe_sparse_dot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mh_bias\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 426\u001b[0;31m                         -(mv.dot(self.W.T)*(mh)).sum(axis=1))     \n\u001b[0m\u001b[1;32m    427\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    428\u001b[0m         Entropy = ( -(mv*np.log(mv)+(1.0-mv)*np.log(1.0-mv)).sum(axis=1)  \n",
      "\u001b[0;32m/Users/charlesmartin14/anaconda/lib/python2.7/site-packages/numpy/core/_methods.pyc\u001b[0m in \u001b[0;36m_sum\u001b[0;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[1;32m     30\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: 'axis' entry is out of bounds"
     ]
    }
   ],
   "source": [
    "rbm._free_energy_TAP(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-10-30T19:34:31.926864",
     "start_time": "2016-10-30T19:34:31.922071"
    },
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "46.0967018276 44.3172430764 90.413944904\n"
     ]
    }
   ],
   "source": [
    "vv = np.dot(v,rbm.v_bias) \n",
    "ff = (np.logaddexp(0, np.dot(v, rbm.W.T)) + rbm.h_bias).sum(axis=0)\n",
    "print vv, ff, vv+ff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# correct - implementation is screwy\n",
    "# fix fix fix"
   ]
  },
  {
   "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
}