## Version 1.000
##
## Code provided by Computational Cognitive Neuroscience Lab
## University of Padova
## ccnl.psy.unipd.it
##
## Permission is granted for anyone to copy, use, modify, or distribute this
## program and accompanying programs and documents for any purpose, provided
## this copyright notice is retained and prominently displayed, along with
## a note saying that the original programs are available from our
## web page.
## The programs and documents are distributed without any warranty, express or
## implied.  As the programs were written for research purposes only, they have
## not been tested to the degree that would be advisable in any important
## application.  All use of these programs is entirely at the user's own risk.
##
## Implementation of a deep belief network on graphic processors (GPUs)


import numpy  as np
import gnumpy as gpu
import scipy.io
import time
gpu.board_id_to_use = 0

# DEEP NETWORK SETUP
# (parameters and final network weights will be saved in structure DN)
DN = {}
DN['layersize'] = np.array([50, 50, 200])
DN['nlayers']   = len(DN['layersize'])
DN['maxepochs'] = 10
DN['batchsize'] = 500
epsilonw_GPU    = 0.1
epsilonvb_GPU   = 0.1
epsilonhb_GPU   = 0.1
weightcost_GPU  = 0.0002
initialmomentum = 0.5
finalmomentum   = 0.9
momentum_GPU    = initialmomentum

dataset = scipy.io.loadmat('MNIST_data_' + str(DN['batchsize']) + '.mat')
batchdata = gpu.garray(dataset['batchdata'])
print '\nUnsupervised pre-training of a deep belief net\n'
x1 = time.strftime('%s')

for layer in range(0, DN['nlayers']):

    # for the first layer, input data are raw images
    # for next layers, input data are preceding hidden activations
    print 'Training layer ', layer+1, '...'
    if layer == 0:
        data_GPU = batchdata
    else:
        data_GPU = batchposhidprobs

    numhid = DN['layersize'][layer]
    numcases, numdims, numbatches = data_GPU.shape
    vishid_GPU       = 0.1 * gpu.randn(numdims, numhid)
    hidbiases_GPU    = gpu.zeros(numhid)
    visbiases_GPU    = gpu.zeros(numdims)
    vishidinc_GPU    = gpu.zeros((numdims, numhid))
    hidbiasinc_GPU   = gpu.zeros(numhid)
    visbiasinc_GPU   = gpu.zeros(numdims)
    batchposhidprobs = gpu.zeros((numcases, numhid, numbatches))

    for epoch in range(DN['maxepochs']):
        for mb in range(numbatches):
            data_mb = data_GPU[:, :, mb]

            #%%%%%%%% START POSITIVE PHASE %%%%%%%%%
            poshidprobs_GPU = (gpu.dot(data_mb, vishid_GPU) + hidbiases_GPU).logistic()
            posprods_GPU    = gpu.dot(data_mb.T, poshidprobs_GPU)
            poshidact_GPU   = poshidprobs_GPU.sum(0)
            posvisact_GPU   = data_GPU[:, :, mb].sum(0)
            #%%%%%%%% END OF POSITIVE PHASE %%%%%%%%%
            poshidstates_GPU = poshidprobs_GPU > poshidprobs_GPU.rand()

            #%%%%%%%% START NEGATIVE PHASE  %%%%%%%%%
            negdata_GPU     = (gpu.dot(poshidstates_GPU, vishid_GPU.T) + visbiases_GPU).logistic()
            neghidprobs_GPU = (gpu.dot(negdata_GPU, vishid_GPU) + hidbiases_GPU).logistic()
            negprods_GPU    = gpu.dot(negdata_GPU.T, neghidprobs_GPU)
            neghidact_GPU   = neghidprobs_GPU.sum(0)
            negvisact_GPU   = negdata_GPU.sum(0)
            #%%%%%%%% END OF NEGATIVE PHASE %%%%%%%%%

            if epoch > 5:
                momentum_GPU = finalmomentum

            #%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%
            vishidinc_GPU  = vishidinc_GPU  * momentum_GPU + ((posprods_GPU - negprods_GPU) / numcases - weightcost_GPU * vishid_GPU) * epsilonw_GPU
            visbiasinc_GPU = visbiasinc_GPU * momentum_GPU + (posvisact_GPU - negvisact_GPU) * (epsilonvb_GPU / numcases)
            hidbiasinc_GPU = hidbiasinc_GPU * momentum_GPU + (poshidact_GPU - neghidact_GPU) * (epsilonhb_GPU / numcases)
            vishid_GPU     = vishid_GPU + vishidinc_GPU
            visbiases_GPU  = visbiases_GPU + visbiasinc_GPU
            hidbiases_GPU  = hidbiases_GPU + hidbiasinc_GPU
            #%%%%%%%% END OF UPDATES %%%%%%%%%

            if epoch == DN['maxepochs']-1:
                batchposhidprobs[:, :, mb] = poshidprobs_GPU

    # save learned weights
    DN['vis_bias' + str(layer)] = visbiases_GPU.as_numpy_array()
    DN['hid_bias' + str(layer)] = hidbiases_GPU.as_numpy_array()
    DN['vishid'   + str(layer)] = vishid_GPU.as_numpy_array()

x2 = time.strftime('%s')
timediff = int(x2) - int(x1)
print '\nElapsed time: ', timediff, ' seconds.'
DN['learningtime'] = timediff
# save final network and parameters
scipy.io.savemat('DN.mat', {'DN': DN})
