- added Rbm

git-svn-id: http://moon:8086/svn/software/trunk/projects/Rbm@560 b431acfa-c32f-4a4a-93f1-934dc6c82436

source/Rbm.cpp

@@ -0,0 +1,228 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+/* 
+ * File:   Rbm.cpp
+ * Author: jens
+ * 
+ * Created on 21. Oktober 2019, 21:28
+ */
+
+#include <streambuf>
+#include "Rbm.hpp"
+#include "noise.h"
+
+Rbm::Rbm(size_t numVisible, size_t numHidden)
+: m_w(numVisible, numHidden)
+, m_bh(1, numHidden)
+, m_bv(1, numVisible)
+{
+	Noise_Init(&m_noise, 0x32727155);
+	uniform(m_w);
+	uniform(m_bh);
+	uniform(m_bv);
+}
+
+Rbm::Rbm(const Rbm& orig)
+{
+}
+
+Rbm::~Rbm()
+{
+}
+
+void Rbm::train(const arma::mat& batch, size_t numEpochs, size_t miniBatchSize, const Params& params, IRbmListener* pListener)
+{
+	size_t i;
+	size_t epoch;
+	size_t gibbs;
+
+	size_t trainingSize = batch.n_rows;
+	size_t trainingSizeRemain = trainingSize;
+	size_t batchRowIndex = 0;
+	
+	double dProgress = 1.0/(numEpochs*(double)trainingSize/std::min(miniBatchSize, trainingSize));
+
+	arma::mat grad_bias_v(arma::zeros(1, m_w.n_rows));
+	arma::mat grad_bias_h(arma::zeros(1, m_w.n_cols));
+	arma::mat grad_weight(arma::zeros(m_w.n_rows, m_w.n_cols));
+	arma::mat momentum_weights = arma::zeros(m_w.n_rows, m_w.n_cols);
+	arma::mat momentum_bias_v(arma::zeros(1, m_w.n_rows));
+	arma::mat momentum_bias_h(arma::zeros(1, m_w.n_cols));
+	arma::mat penalty_weights = arma::zeros(m_w.n_rows, m_w.n_cols);
+	double L1 = 0;
+	double L2 = 0;
+	
+	double progress = 0;
+	while (trainingSizeRemain)
+	{
+		std::cout << "trainingSizeRemain: " << trainingSizeRemain << std::endl;
+		size_t toSlice = std::min(miniBatchSize, trainingSizeRemain);
+		arma::mat miniBatch = batch.rows(batchRowIndex, batchRowIndex+toSlice-1);
+		trainingSizeRemain -= toSlice;
+		batchRowIndex += toSlice;
+		size_t miniBatchSizeActual = miniBatch.n_rows;
+		double learning_rate = params.m_learningRate/std::min(miniBatchSizeActual, trainingSize);
+		double weight_decay = params.m_weightDecay/std::min(miniBatchSizeActual, trainingSize);
+
+		arma::mat vis_state(miniBatchSizeActual, m_w.n_rows);
+		arma::mat vis_probs(miniBatchSizeActual, m_w.n_rows);
+		arma::mat hid_state(miniBatchSizeActual, m_w.n_cols);
+		arma::mat hid_probs(miniBatchSizeActual, m_w.n_cols);
+
+		for (epoch=0; epoch < numEpochs; epoch++)
+		{
+			if (pListener)
+			{
+				if(!pListener->onProgress())
+				{
+					break;
+				}
+			}
+			
+			// Create hidden layer base on training data
+			if (params.m_doSampleBatch)
+			{
+				// When the hidden units are being driven by data, always use stochastic binary states
+				vis_state = sample(miniBatch);
+			}
+			else
+			{
+				vis_state = miniBatch;
+			}
+			
+			hid_state = vis_state * m_w + arma::repmat(m_bh, miniBatchSizeActual, 1);
+			hid_probs = probsLogistic(hid_state);
+
+			// Sample hidden
+			if (params.m_doRaoBlackwell)
+			{
+				hid_state = hid_probs;
+			}
+			else
+			{
+				hid_state = sample(hid_probs);
+			}
+
+			// Update weights (positive phase)
+			grad_weight = vis_state.t() * hid_state;
+			grad_bias_v = sum(vis_state, 0);
+			grad_bias_h = sum(hid_state, 0);
+
+			for (gibbs=0; gibbs < params.m_numGibbs; gibbs++)
+			{
+
+				// Create hidden representation given v
+				hid_state = sample(hid_probs);
+
+				// Create visible reconstruction (a fantasy...) given hid
+				vis_state = hid_state * m_w.t() + arma::repmat(m_bv, miniBatchSizeActual, 1);
+				vis_probs = probsLogistic(vis_state);
+				if (params.m_doSampleVisible)
+				{
+					vis_state = sample(vis_probs);
+				}
+				else
+				{
+					vis_state = vis_probs;
+				}
+
+				// Create hidden representation given v
+				hid_state = vis_state * m_w + repmat(m_bh, miniBatchSizeActual, 1);
+				hid_probs = probsLogistic(hid_state);
+			
+			}
+
+			// Update weights (negative phase)
+			grad_bias_v -= sum(vis_probs, 0);
+			grad_bias_h -= sum(hid_probs, 0);
+			grad_weight -= vis_probs.t() * hid_probs;
+
+			for (int i=0; i < m_w.n_rows; i++)
+			{
+				for (int j=0; j < m_w.n_cols; j++)
+				{
+					if (m_w(i,j) >= 0)
+					{
+						penalty_weights(i,j) = weight_decay;
+					}
+					else
+					{
+						penalty_weights(i,j) = -weight_decay;
+					}
+				}
+			}
+			
+			L1 = accu(abs(m_w));
+			L2 = accu(m_w % m_w);
+			momentum_bias_v = params.m_momentum*momentum_bias_v + grad_bias_v;
+			momentum_bias_h = params.m_momentum*momentum_bias_h + grad_bias_h;
+			momentum_weights = params.m_momentum*momentum_weights + grad_weight - L2*penalty_weights;
+			
+			m_bv += learning_rate*momentum_bias_v;
+			m_bh += learning_rate*momentum_bias_h;
+			m_w += learning_rate*momentum_weights;
+			
+			progress += dProgress;
+
+		} // Number of epochs
+			
+		arma::mat diffErr = miniBatch - vis_probs;
+		arma::mat diffErr_squared = diffErr % diffErr;
+		double err = accu(diffErr_squared);
+		std::cout << "error (per mini batch) = " << err << std::endl;
+		std::cout << "L1                     = " << L1 << std::endl;
+		std::cout << "L2                     = " << L2 << std::endl;
+	} // number of mini batches
+}
+
+arma::mat Rbm::probsLogistic(const arma::mat &src)
+{
+	return 1 / (1 + (arma::exp(-src)));
+}
+
+arma::mat Rbm::sample(const arma::mat &src)
+{
+	arma::mat dst = src;
+	for (size_t i=0; i < src.n_rows; i++)
+	{
+		for (size_t j=0; j < src.n_cols; j++)
+		{
+			dst(i, j) = src(i, j) >= Noise_Uniform(&m_noise);
+		}
+	}
+	return dst;
+}
+
+arma::mat Rbm::toHidden(const arma::mat &visible)
+{
+	arma::mat h = visible.t() * m_w  + m_bh;
+	return h;
+}
+
+arma::mat Rbm::toVisible(const arma::mat &hidden)
+{
+	arma::mat v = hidden * m_w.t() + m_bv;
+	return v;
+}
+
+arma::mat Rbm::uniform(size_t numRows, size_t numCols, double mu, double stdDev)
+{
+	arma::mat dst = arma::zeros<arma::mat>(numRows, numCols);
+	uniform(dst);
+	return dst;
+}
+
+void Rbm::uniform(arma::mat& srcDst, double mu, double stdDev)
+{
+	for (size_t i=0; i < srcDst.n_rows; i++)
+	{
+		for (size_t j=0; j < srcDst.n_cols; j++)
+		{
+			srcDst(i, j) = stdDev*Noise_Uniform(&m_noise) + mu;
+		}
+	}
+}

source/Rbm.hpp

@@ -0,0 +1,77 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+/* 
+ * File:   Rbm.hpp
+ * Author: jens
+ *
+ * Created on 21. Oktober 2019, 21:28
+ */
+
+#ifndef RBM_HPP
+#define RBM_HPP
+
+#include <armadillo>
+#include "noise.h"
+class IRbmListener
+{
+public:
+	
+	virtual ~IRbmListener()
+	{
+
+	}
+
+	bool onProgress()
+	{
+		return true;
+	}
+};
+class Rbm
+{
+public:
+	struct Params
+	{
+		Params()
+		: m_weightDecay(0.01)
+		, m_learningRate(0.1)
+		, m_momentum(0.5)
+		, m_doRaoBlackwell(true)
+		, m_doSampleVisible(false)
+		, m_doSampleBatch(false)
+		, m_numGibbs(1)
+		{
+		}
+		
+		double m_weightDecay;
+		double m_learningRate;
+		double m_momentum;
+		bool m_doRaoBlackwell;
+		bool m_doSampleVisible;
+		bool m_doSampleBatch;
+		size_t m_numGibbs;
+	};
+
+	Rbm(size_t numHidden, size_t numVisible);
+	Rbm(const Rbm& orig);
+	virtual ~Rbm();
+
+	void train(arma::mat const &batch, size_t numEpochs, size_t sizeMiniBatch, Params const &params, IRbmListener *pListener);
+	arma::mat sample(arma::mat const &src);
+	static arma::mat probsLogistic(arma::mat const &src);
+	arma::mat toHidden(const arma::mat &v);
+	arma::mat toVisible(const arma::mat &h);
+	arma::mat uniform(size_t numRows, size_t numCols, double mu=0.0, double stdDev=1.0);
+	void uniform(arma::mat &srcDst, double mu=0.0, double stdDev=1.0);
+private:
+	noise_gen_t m_noise;
+	arma::mat m_w;
+	arma::mat m_bh;
+	arma::mat m_bv;
+};
+
+#endif /* RBM_HPP */
+

source/main.cpp

@@ -1,6 +1,7 @@
 #include <cstdio>
 #include <cmath>
 #include <armadillo>
+#include "Rbm.hpp"
 
 int main()
 {
@@ -23,12 +24,20 @@ int main()
 	Pos.print("New position of the particle:");   //    ^
 												//    x (1,0)
 												//    | 
-	arma::mat Z = arma::eye<arma::mat>(4,4);
+	arma::mat Z = arma::eye<arma::mat>(4000,4000);
-	Z.print("Z:");
+//	Z.print("Z:");
 	printf("Z.n_rows = %d\n", (int)Z.n_rows);
 	printf("Z.n_cols = %d\n", (int)Z.n_cols);
 	printf("Z.n_elem = %d\n", (int)Z.n_elem);
 	//    +------>
 
+	Rbm::Params params;
+	Rbm rbm(28*28, 64);
+	arma::mat batch = arma::randu(1000, 28*28);
+	rbm.train(batch, 1000, 100, params, nullptr);
+	
+//	arma::mat v = arma::randu(28*28, 1);
+//	arma::mat h = rbm.toHidden(v);
+//	arma::mat r = rbm.toVisible(h);
 	return 0;
 }

source/noise.c

@@ -0,0 +1,80 @@
+// --------------------------------------------------------------
+// --------------------------------------------------------------
+
+#include <string.h>
+#include <stdlib.h>
+#include <time.h>
+#include <math.h>
+
+#include "noise.h"
+
+// --------------------------------------------------------------
+// internal funcs
+// --------------------------------------------------------------
+#define PM_IA 16807
+#define PM_IM 2147483647
+#define PM_AM (1.0/PM_IM)
+#define PM_IQ 127773
+#define PM_IR 2836
+#define PM_MASK 123459876
+
+//	'Minimal' random number generator of Park and Miller. Returns a uniform random deviate
+//	between 0.0 and 1.0. Set or reset idum to any integer value (except the unlikely value MASK)
+//	to initialize the sequence; idum must not be altered between calls for successive deviates in
+//	a sequence.
+double ran0(long *idum)
+{
+	long k;
+	double ans;
+
+	*idum ^= PM_MASK;				// XORing with MASK allows use of zero and other
+	k=(*idum)/PM_IQ;				// simple bit patterns for idum.
+	*idum=PM_IA*(*idum-k*PM_IQ)-PM_IR *k;		// Compute idum=(IA*idum) % IM without over-
+							// flows by Schrage’s method.
+	if (*idum < 0)
+		*idum += PM_IM;
+
+	ans=PM_AM*(*idum);				// Convert idumto a floating result.
+	*idum ^= PM_MASK;				// Unmask before return.
+
+	return ans;
+}
+
+// --------------------------------------------------------------
+// Exported functions
+// --------------------------------------------------------------
+void Noise_Init(noise_gen_t *pObj, long seed)
+{
+
+	pObj->state = seed ^ clock();
+}
+
+void Noise_Free(noise_gen_t *pObj)
+{
+}
+
+// 0 .. 1
+double Noise_Uniform(noise_gen_t *pObj)
+{
+	return ran0(&pObj->state);
+}
+
+double Noise_Gaussian(noise_gen_t *pObj)
+{
+	double U1, U2, V1, V2, S, Y;
+
+	do
+	{
+		U1 = Noise_Uniform(pObj);		         // U1=[0,1]
+		U2 = Noise_Uniform(pObj);                // U2=[0,1]
+		V1 = 2 * U1 - 1;          // V1=[-1,1]
+		V2 = 2 * U2 - 1;          // V2=[-1,1]
+		S = V1 * V1 + V2 * V2;
+
+	} while (S >= 1);
+
+	// X = sqrt(-2 * log(S) / S) * V1;
+	Y = sqrt(-2 * log(S) / S) * V2;
+
+	return Y;
+}

source/noise.h

@@ -0,0 +1,42 @@
+/*
+  ==============================================================================
+
+    noise.h
+    Created: 21 Sep 2014 1:55:07pm
+    Author:  jens
+
+  ==============================================================================
+*/
+// --------------------------------------------------------------
+#ifndef _NOISE_H_
+#define	_NOISE_H_
+
+#define __func__ ""
+// --------------------------------------------------------------
+// Types
+// --------------------------------------------------------------
+typedef struct _snoise_gen_t
+{
+	long state;
+} noise_gen_t;
+
+// --------------------------------------------------------------
+#if defined(__cplusplus)
+extern "C" {
+#endif
+
+// --------------------------------------------------------------
+// Exported functions
+// --------------------------------------------------------------
+void Noise_Init(noise_gen_t *pObj, long seed);
+void Noise_Free(noise_gen_t *pObj);
+double Noise_Uniform(noise_gen_t *pObj);
+double Noise_Gaussian(noise_gen_t *pObj);
+
+#if defined(__cplusplus)
+}
+
+#endif
+// --------------------------------------------------------------
+
+#endif // _NOISE_H_