RBMPaper-Rcode/usageExample.r at main · ph09/RBMPaper-Rcode · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
rm(list=ls())
source('rbm.r')
source('readData.r')
source('obtainHiddenRep.r')
source('forward.r')
source('forwardPass.r')
source('sigmoid.r')
source('softThresholding.r')
graphics.off()
##read data
library(R.matlab)
library(lattice)

datasetName = 'dataset1.mat' #path to dataset file.
#the dataset needs to be a .mat file, with a binary matrix f of size d x n
#and (optionally, a binary label vector y).
data = readData(datasetName)

## setup
rbmInput <- list()
rbmOutput <- list()
monitorInput <- list()
rbmInput$restart=1;

#the following are configurable hyperparameters for RBM

rbmInput$reg_type = 'l2'
rbmInput$weightPenalty = 1e-2 #\ell_2 weight penalty
weightPenaltyOrg = rbmInput$weightPenalty
rbmInput$epsilonw      = 5e-2 #5e-2 # Learning rate for weights
rbmInput$epsilonvb     = 5e-2 #5e-2 # Learning rate for biases of visible units
rbmInput$epsilonhb     = 5e-2 #5e-2 # Learning rate for biases of hidden units
rbmInput$CD=10 #number of contrastive divergence iterations
rbmInput$initialmomentum  = 0
rbmInput$finalmomentum    = 0.9
rbmInput$maxEpoch = 150
rbmInput$decayLrAfter = 120
rbmInput$decayMomentumAfter = 90 #when to switch from initial to final momentum
rbmInput$iIncreaseCD = 0
# monitor free energy and likelihood change (on validation set) with time
rbmInput$iMonitor = 1


## train
sizes = list()
rbmInput$data = data
rbmInput$numhid = ncol(data$allDataTable)
stack = vector('list', 1)
layerCounter = 1
addLayers = 1
while (addLayers){
  # train RBM
  rbmInput$weightPenalty = weightPenaltyOrg
  rbmOutput = rbm(rbmInput)
  # collect params
  stack[[layerCounter]]$vishid = rbmOutput$vishid
  stack[[layerCounter]]$hidbiases = rbmOutput$hidbiases
  stack[[layerCounter]]$visbiases = rbmOutput$visbiases

  # SVD to determine number of hidden nodes
  tmp = svd (stack[[layerCounter]]$vishid)
  U = tmp$u
  D = diag(tmp$d)
  V = tmp$v
  numhid = min(which((cumsum(D)/sum(D))>0.95))
  message(sprintf ('need %1.0f hidden units\n', numhid))
  print('paused, press Enter key to continue')
  scan(quiet=TRUE)

  # Re-train RBM
  sizes = c(sizes, numhid)
  rbmInput$numhid = numhid
  rbmInput$weightPenalty = 0 #rbmInput.weightPenalty/10;
  rbmOutput = rbm(rbmInput)
  # collect params
  stack[[layerCounter]]$vishid = rbmOutput$vishid
  stack[[layerCounter]]$hidbiases = rbmOutput$hidbiases
  stack[[layerCounter]]$visbiases = rbmOutput$visbiases
  v=c('weight matrix of RBM ', as.character(layerCounter))
  v = paste(v, collapse='')
  #levelplot(stack[[layerCounter]]$vishid, col.regions=terrain.colors(100),
            #scales=list(x=list(0:((ncol(stack[[layerCounter]]$vishid)>5)+1):ncol(stack[[layerCounter]]$vishid)),
            #y=list(0:((nrow(stack[[layerCounter]]$vishid)>5)+1):nrow(stack[[layerCounter]]$vishid))),
            #main=v, xlab='hidden units', ylab='visible units')
  # setup for next RBM
  rbmInput$data = obtainHiddenRep(rbmInput, rbmOutput)

  # stopping criterion
  if (numhid ==1){
    addLayers = 0
  }
  layerCounter = layerCounter + 1
}

numLayers = ncol(stack)
message(sprintf('trained a deep net with %1.0f layers, of sizes:\n', numLayers))
print(sizes)
## obtain posterior probabilities
#deterministic
mode = 'deterministic'
posteriorProbsDet = forward (stack, data$allDataTable, mode)


# stochastic
mode = 'stochastic'
nit = 100
posteriorProbsStoch = forward(stack, data$allDataTable, mode, nit)

## predict labels
labels = data$labels

# deterministic mode:
predictedLabels = round(posteriorProbsDet)
# check if predictedLables need to be flipped
m = mean(predictedLabels == data$allDataTable[, 1])
if (m<0.5){
  predictedLabels = 1 - predictedLabels
}
acc = mean(labels==predictedLabels)
inds1 = labels[labels==1]
inds0 = labels[labels==0]
sensitivity = mean(predictedLabels[inds1])
specificity = 1-mean(predictedLabels[inds0])

balAcc_rbmDet = (sensitivity + specificity)/2
print('Deterministic mode:')
message(sprintf(1,'sensitivity: %0.3f%%\n',100*sensitivity))
message(sprintf(1,'specificity: %0.3f%%\n',100*specificity))
message(sprintf(1,'accuracy: %0.3f%%\n',100*acc))
message(sprintf(1,'balanced accuracy: %0.3f%%\n',100*balAcc_rbmDet))

#stochastic mode:
predictedLabels = round(posteriorProbsStoch)
#check if predictedLables need to be flipped
m = mean(predictedLabels == data$allDataTable[, 1])
if (m<0.5){
  predictedLabels = 1-predictedLabels
}
acc = mean(labels==predictedLabels)
inds1 = labels[labels==1]
inds0 = labels[labels==0]
sensitivity = mean(predictedLabels[inds1])
specificity = 1-mean(predictedLabels[inds0])

balAcc_rbmStoch = (sensitivity + specificity)/2
print('Stochastic mode:')
message(sprintf(1,'sensitivity: %0.3f%%\n',100*sensitivity))
message(sprintf(1,'specificity: %0.3f%%\n',100*specificity))
message(sprintf(1,'accuracy: %0.3f%%\n',100*acc))
message(sprintf(1,'balanced accuracy: %0.3f%%\n',100*balAcc_rbmStoch))