<commit_before>"""
This tutorial introduces restricted boltzmann machines (RBM) using Theano.

Description of RBM : TODO
"""


import numpy
import theano
import theano.tensor as T
import time
import gzip
import cPickle

from theano.tensor.shared_randomstreams import RandomStreams

from theano.sandbox.scan import scan


class RBM():
    """Restricted Boltzmann Machine (RBM)
    """
    def __init__(self, input=None, n_visible=784, n_hidden=500, \
        batch_size = 20, shared_W = None, shared_b = None):
        """ 
        RBM constructor. Defines the parameters of the model along with
        basic operations for inferring hidden from visible (and vice-versa), 
        as well as for performing CD updates.

        :param input: None for standalone RBMs or symbolic variable if RBM is
        part of a larger graph.

        :param n_visible: number of visible units

        :param n_hidden: number of hidden units

        :param batch_size: [mini]batch size used in the SGD updates

        :param shared_W: None for standalone RBMs or symbolic variable to a
        shared weight matrix in case RBM is part of a DBN network; in a DBN,
        the weights are shared between RBMs and layers of a MLP

        :param shared_b: None for standalone RBMs or symbolic variable to a
        shared bias vector in case RBM is part of a DBN network
        """

        self.n_visible = n_visible
        self.n_hidden  = n_hidden
        
        # setup theano random number generator
        self.theano_rng = RandomStreams()
        numpy_rng  = numpy.random.RandomState()
        
        # initial values for weights and biases
        if shared_W and shared_b :
            self.W     = shared_W
            self.hbias = shared_b
        else:
            # W is initialized with `initial_W` which is uniformely sampled
            # from -6./sqrt(n_visible+n_hidden) and 6./sqrt(n_hidden+n_visible)
            # the output of uniform if converted using asarray to dtype 
            # theano.config.floatX so that the code is runable on GPU
            initial_W = numpy.asarray( numpy.random.uniform( \
                  low = -numpy.sqrt(6./(n_hidden+n_visible)), \
                high = numpy.sqrt(6./(n_hidden+n_visible)), \
                size = (n_visible, n_hidden)), dtype = theano.config.floatX)
            # initial value of the hidden units bias
            initial_hbias   = numpy.zeros(n_hidden)

            # theano shared variables for weights and biases
            self.W     = theano.shared(value = initial_W    , name = 'W')
            self.hbias = theano.shared(value = initial_hbias, name = 'hbias')


        # initial value of the visible units bias 
        initial_vbias  = numpy.zeros(n_visible)
        self.vbias = theano.shared(value = initial_vbias, name = 'vbias')

        # initialize input layer for standalone RBM or layer0 of DBN
        self.input = input if input else T.dmatrix('input')

        # Create the shared variables that will store the negative phase
        # values; these values are usable only in case of PCD
        self.n_vis = theano.shared(value = numpy.zeros((batch_size,n_visible)))
        self.n_hid = theano.shared(value = numpy.zeros((batch_size,n_hidden)))
        self.params     = [self.W, self.vbias, self.hbias]
        self.batch_size = batch_size
 
    def CD_k(self, n_Gibbs_steps = 1, persistent = False):
        '''
        :param n_Gibbs_steps: number of Gibbs steps to do when computing the 
        gradient with CD
        '''

        # define the graph for positive phase
        p_hid_activation = T.dot(self.input, self.W) + self.hbias
        p_hid            = T.nnet.sigmoid(p_hid_activation)
        p_hid_sample     = self.theano_rng.binomial(T.shape(p_hid), 1, p_hid)*1.0
        
        # for negative phase we need to implement k Gibbs steps; for this we 
        # will use the scan op (see theano documentation about it)
    
        # Create a function that builds the graph for one step of Gibbs
        def oneGibbsStep( vis_km1, hid_km1):
           # For clarity we used the following naming conventions : 
           #     var_kmx -> the value of variable ``var`` at step k - x
           #     var_k   -> the value of variable ``var`` at step k
           #     var_kpx -> the value of variable ``var`` at step k + x

           # if you want persistent CD, we have to use the last negative 
           # value generated ( which is stored in ``self.n_hid``
           if persistent : 
                hid_km1 = self.n_hid

           # sample this hidden layer
           hid_sample_km1   = self.theano_rng.binomial(T.shape(hid_km1),1,hid_km1)*1.0
           # compute visible layer values
           vis_activation_k = T.dot(hid_sample_km1, self.W.T) + self.vbias
           vis_k            = T.nnet.sigmoid(vis_activation_k)
           # sample the visible
           vis_sample_k     = self.theano_rng.binomial(T.shape(vis_k),1,vis_k)*1.0
           # compute new hidden layer values
           hid_activation_k = T.dot(vis_sample_k, self.W) + self.hbias
           hid_k            = T.nnet.sigmoid(hid_activation_k)

           # return a list of outputs, plus a dictionary of updates 
           return ([vis_k, hid_k],{ self.n_vis : vis_k, self.n_hid : hid_k})
        
        # to compute the negative phase perform k Gibbs step; for this we 
        # use the scan op, that implements a loop

        # keep_outputs tells scan that we do not care about intermediate values
        # of n_vis and n_hid, and that it should only return the last one
        n_vis, n_hid = scan(oneGibbsStep,[],[self.input, p_hid],\
               [], n_steps = n_Gibbs_steps, keep_outputs = {0:False, 1:False} ) 

        g_vbias = T.mean( self.input - n_vis , axis = 0)
        g_hbias = T.mean( p_hid      - n_hid , axis = 0)
        
        g_W = T.dot(p_hid.T, self.input )/ self.batch_size - \
              T.dot(n_hid.T, n_vis      )/ self.batch_size

        gparams = [g_W.T, g_vbias, g_hbias]
        # define dictionary of stochastic gradient update equations
        cost = T.mean(abs(self.input - n_vis))
        return (gparams, cost)

       


def sgd_optimization_mnist( learning_rate=0.1, training_epochs = 20, \
                            dataset='mnist.pkl.gz'):
    """
    Demonstrate stochastic gradient descent optimization for a multilayer 
    perceptron

    This is demonstrated on MNIST.

    :param learning_rate: learning rate used in the finetune stage 
    (factor for the stochastic gradient)

    :param pretraining_epochs: number of epoch to do pretraining

    :param pretrain_lr: learning rate to be used during pre-training

    :param n_iter: maximal number of iterations ot run the optimizer <commit_msg>comments on rbm.py<commit_after>"""
This tutorial introduces restricted boltzmann machines (RBM) using Theano.

Description of RBM : TODO
"""


import numpy
import theano
import theano.tensor as T
import time
import gzip
import cPickle

from theano.tensor.shared_randomstreams import RandomStreams

from theano.sandbox.scan import scan


class RBM():
    """Restricted Boltzmann Machine (RBM)
    """
    def __init__(self, input=None, n_visible=784, n_hidden=500, \
        batch_size = 20, shared_W = None, shared_b = None):
        """ 
        RBM constructor. Defines the parameters of the model along with
        basic operations for inferring hidden from visible (and vice-versa), 
        as well as for performing CD updates.

        :param input: None for standalone RBMs or symbolic variable if RBM is
        part of a larger graph.

        :param n_visible: number of visible units

        :param n_hidden: number of hidden units

        :param batch_size: [mini]batch size used in the SGD updates

        :param shared_W: None for standalone RBMs or symbolic variable to a
        shared weight matrix in case RBM is part of a DBN network; in a DBN,
        the weights are shared between RBMs and layers of a MLP

        :param shared_b: None for standalone RBMs or symbolic variable to a
        shared bias vector in case RBM is part of a DBN network
        """

        self.n_visible = n_visible
        self.n_hidden  = n_hidden
        
        # setup theano random number generator
        self.theano_rng = RandomStreams()
        numpy_rng  = numpy.random.RandomState()
        
        # initial values for weights and biases
        if shared_W and shared_b :
            self.W     = shared_W
            self.hbias = shared_b
        else:
            # W is initialized with `initial_W` which is uniformely sampled
            # from -6./sqrt(n_visible+n_hidden) and 6./sqrt(n_hidden+n_visible)
            # the output of uniform if converted using asarray to dtype 
            # theano.config.floatX so that the code is runable on GPU
            initial_W = numpy.asarray( numpy.random.uniform( \
                  low = -numpy.sqrt(6./(n_hidden+n_visible)), \
                high = numpy.sqrt(6./(n_hidden+n_visible)), \
                size = (n_visible, n_hidden)), dtype = theano.config.floatX)
            # initial value of the hidden units bias
            initial_hbias   = numpy.zeros(n_hidden)

            # theano shared variables for weights and biases
            self.W     = theano.shared(value = initial_W    , name = 'W')
            self.hbias = theano.shared(value = initial_hbias, name = 'hbias')


        # initial value of the visible units bias 
        initial_vbias  = numpy.zeros(n_visible)
        self.vbias = theano.shared(value = initial_vbias, name = 'vbias')

        # initialize input layer for standalone RBM or layer0 of DBN
        self.input = input if input else T.dmatrix('input')

        # Create the shared variables that will store the negative phase
        # values; these values are usable only in case of PCD
        self.n_vis = theano.shared(value = numpy.zeros((batch_size,n_visible)))
        self.n_hid = theano.shared(value = numpy.zeros((batch_size,n_hidden)))
        self.params     = [self.W, self.vbias, self.hbias]
        self.batch_size = batch_size
 
    def CD_k(self, n_Gibbs_steps = 1, persistent = False):
        '''
        :param n_Gibbs_steps: number of Gibbs steps to do when computing the 
        gradient with CD

        *** Document what does this function return?

        *** How about a more flexible way to make this function do PCD: an optional
        visible_0=None parameter to this function where None means self.input.  This variable
        is used to initialize the negative chain.
        '''

        # define the graph for positive phase
        p_hid_activation = T.dot(self.input, self.W) + self.hbias
        p_hid            = T.nnet.sigmoid(p_hid_activation)
        p_hid_sample     = self.theano_rng.binomial(T.shape(p_hid), 1, p_hid)*1.0
        
        # for negative phase we need to implement k Gibbs steps; for this we 
        # will use the scan op (see theano documentation about it)
    
        # Create a function that builds the graph for one step of Gibbs
        def oneGibbsStep( vis_km1, hid_km1):
           # For clarity we used the following naming conventions : 
           #     var_kmx -> the value of variable ``var`` at step k - x
           #     var_k   -> the value of variable ``var`` at step k
           #     var_kpx -> the value of variable ``var`` at step k + x

           # if you want persistent CD, we have to use the last negative 
           # value generated ( which is stored in ``self.n_hid``

           # *** 
           if persistent : 
                hid_km1 = self.n_hid

           # sample this hidden layer
           hid_sample_km1   = self.theano_rng.binomial(T.shape(hid_km1),1,hid_km1)*1.0
           # compute visible layer values
           vis_activation_k = T.dot(hid_sample_km1, self.W.T) + self.vbias
           vis_k            = T.nnet.sigmoid(vis_activation_k)
           # sample the visible
           vis_sample_k     = self.theano_rng.binomial(T.shape(vis_k),1,vis_k)*1.0
           # compute new hidden layer values
           hid_activation_k = T.dot(vis_sample_k, self.W) + self.hbias
           hid_k            = T.nnet.sigmoid(hid_activation_k)

           # return a list of outputs, plus a dictionary of updates 

           # *** Could we modify scan to accept a friendlier encoding of this information,
           # like: 
           #    dict( outputs=[vis_k, hid_k], updates={...})
           return ([vis_k, hid_k],{ self.n_vis : vis_k, self.n_hid : hid_k})
        
        # to compute the negative phase perform k Gibbs step; for this we 
        # use the scan op, that implements a loop


        # *** Could we use a different variable name here?  n_something usually means the
        # number of somethings.  Like the number of visible or hidden units.

        # keep_outputs tells scan that we do not care about intermediate values
        # of n_vis and n_hid, and that it should only return the last one
        n_vis, n_hid = scan(oneGibbsStep,[],[self.input, p_hid],\
               [], n_steps = n_Gibbs_steps, keep_outputs = {0:False, 1:False} ) 

        g_vbias = T.mean( self.input - n_vis , axis = 0)
        g_hbias = T.mean( p_hid      - n_hid , axis = 0)

        # ***Why are we using mean for the biases but a dot()/size formula for the weights?
        #    It's a minor point, but we're confusing two kinds of terminology.
        #    Better would be using mean & covariance (I think we have a cov() op...)
        #    -or- sum()/batchsize and then dot() / batchsize
        
        g_W = T.dot(p_hid.T, self.input )/ self.batch_size - \
              T.dot(n_hid.T, n_vis      )/ self.batch_size

        gparams = [g_W.T, g_vbias, g_hbias]
        # define dictionary of stochastic gradient update equations

        # *** It is misleading to say that it returns a cost, since usually a cost is the thing
        # we are minimizing.

        cost = T.mean(abs(self.input - n_vis))
        return (gparams, cost)

       

# *** rename this function
def sgd_optimization_mnist( learning_rate=0.1, training_epochs = 20, \
                            dataset='mnist.pkl.gz'):
    """
    Demonstrate ***

    This is demonstrated on MNIST.

    :param learning_rate: learning rate used in the finetune stage 
    (factor for the stochastic gradient)

    :param pretraining_epochs: number of epoch to do pretraining

    :param pretrain_lr: learning rate to be used during pre-training

    :param n_iter: maximal number of iterations ot run the optimizer 
{'fraction_non_alphanumeric': 0.07650610676547276, 'fraction_numerical': 0.005283381364073006, 'mean_word_length': 2.183009390696659, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4061512', 'n_tokens_mistral': 4407, 'n_tokens_neox': 4219, 'n_words': 1482}