1 import numpy 
  2 from scipy import sparse as S
  3 from matplotlib import pyplot as plt 
  4 from scipy.sparse.csr import csr_matrix 
  5 import pandas 
  6 
  7 def normalize(x):
  8     V = x.copy()
  9     V -= x.min(axis=1).reshape(x.shape[0],1)
 10     V /= V.max(axis=1).reshape(x.shape[0],1)
 11     return V
 12     
 13 def sigmoid(x):
 14     #return x*(x > 0)
 15     #return numpy.tanh(x)
 16     return 1.0/(1+numpy.exp(-x)) 
 17 
 18 class RBM():
 19     def __init__(self, n_visible=None, n_hidden=None, W=None, learning_rate = 0.1, weight_decay=1,cd_steps=1,momentum=0.5):
 20         if W == None:
 21             self.W =  numpy.random.uniform(-.1,0.1,(n_visible,  n_hidden)) / numpy.sqrt(n_visible + n_hidden)
 22             self.W = numpy.insert(self.W, 0, 0, axis = 1)
 23             self.W = numpy.insert(self.W, 0, 0, axis = 0)
 24         else:
 25             self.W=W 
 26         self.learning_rate = learning_rate 
 27         self.momentum = momentum
 28         self.last_change = 0
 29         self.last_update = 0
 30         self.cd_steps = cd_steps
 31         self.epoch = 0 
 32         self.weight_decay = weight_decay  
 33         self.Errors = []
 34          
 35             
 36     def fit(self, Input, max_epochs = 1, batch_size=100):  
 37         if isinstance(Input, S.csr_matrix):
 38             bias = S.csr_matrix(numpy.ones((Input.shape[0], 1))) 
 39             csr = S.hstack([bias, Input]).tocsr()
 40         else:
 41             csr = numpy.insert(Input, 0, 1, 1)
 42         for epoch in range(max_epochs): 
 43             idx = numpy.arange(csr.shape[0])
 44             numpy.random.shuffle(idx)
 45             idx = idx[:batch_size]  
 46                    
 47             self.V_state = csr[idx] 
 48             self.H_state = self.activate(self.V_state)
 49             pos_associations = self.V_state.T.dot(self.H_state) 
 50   
 51             for i in range(self.cd_steps):
 52               self.V_state = self.sample(self.H_state)  
 53               self.H_state = self.activate(self.V_state)
 54               
 55             neg_associations = self.V_state.T.dot(self.H_state) 
 56             self.V_state = self.sample(self.H_state) 
 57             
 58             # Update weights. 
 59             w_update = self.learning_rate * ((pos_associations - neg_associations) / batch_size) 
 60             total_change = numpy.sum(numpy.abs(w_update)) 
 61             self.W += self.momentum * self.last_change  + w_update
 62             self.W *= self.weight_decay 
 63             
 64             self.last_change = w_update
 65             
 66             RMSE = numpy.mean((csr[idx] - self.V_state)**2)**0.5
 67             self.Errors.append(RMSE)
 68             self.epoch += 1
 69             print("Epoch %s: RMSE = %s; ||W||: %6.1f; Sum Update: %f" % (self.epoch, RMSE, numpy.sum(numpy.abs(self.W)), total_change))  
 70         return self 
 71         
 72     def learning_curve(self):
 73         plt.ion()
 74         #plt.figure()
 75         plt.show()
 76         E = numpy.array(self.Errors)
 77         plt.plot(pandas.rolling_mean(E, 50)[50:])  
 78      
 79     def activate(self, X):
 80         if X.shape[1] != self.W.shape[0]:
 81             if isinstance(X, S.csr_matrix):
 82                 bias = S.csr_matrix(numpy.ones((X.shape[0], 1))) 
 83                 csr = S.hstack([bias, X]).tocsr()
 84             else:
 85                 csr = numpy.insert(X, 0, 1, 1) 
 86         else:
 87             csr = X
 88         p = sigmoid(csr.dot(self.W)) 
 89         p[:,0]  = 1.0 
 90         return p  
 91         
 92     def sample(self, H, addBias=True): 
 93         if H.shape[1] == self.W.shape[0]:
 94             if isinstance(H, S.csr_matrix):
 95                 bias = S.csr_matrix(numpy.ones((H.shape[0], 1))) 
 96                 csr = S.hstack([bias, H]).tocsr()
 97             else:
 98                 csr = numpy.insert(H, 0, 1, 1)
 99         else:
100             csr = H
101         p = sigmoid(csr.dot(self.W.T)) 
102         p[:,0] = 1
103         return p
104       
105 if __name__=="__main__":
106     data = numpy.random.uniform(0,1,(100,10))
107     rbm = RBM(10,15)
108     rbm.fit(data,1000)
109     rbm.learning_curve()