From 6ad74f4ee7042f53af81b5c4c3758050c6dab825 Mon Sep 17 00:00:00 2001 From: Jens Ahrensfeld Date: Mon, 21 Oct 2019 21:28:26 +0000 Subject: [PATCH] - added Rbm git-svn-id: http://moon:8086/svn/software/trunk/projects/Rbm@560 b431acfa-c32f-4a4a-93f1-934dc6c82436 --- source/Rbm.cpp | 228 ++++++++++++++++++++++++++++++++++++++++++++++++ source/Rbm.hpp | 77 ++++++++++++++++ source/main.cpp | 13 ++- source/noise.c | 80 +++++++++++++++++ source/noise.h | 42 +++++++++ 5 files changed, 438 insertions(+), 2 deletions(-) create mode 100644 source/Rbm.cpp create mode 100644 source/Rbm.hpp create mode 100644 source/noise.c create mode 100644 source/noise.h diff --git a/source/Rbm.cpp b/source/Rbm.cpp new file mode 100644 index 0000000..3f8e7fd --- /dev/null +++ b/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 +#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(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; + } + } +} diff --git a/source/Rbm.hpp b/source/Rbm.hpp new file mode 100644 index 0000000..ce33454 --- /dev/null +++ b/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 +#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 ¶ms, 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 */ + diff --git a/source/main.cpp b/source/main.cpp index e61efef..423e15c 100644 --- a/source/main.cpp +++ b/source/main.cpp @@ -1,6 +1,7 @@ #include #include #include +#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(4,4); - Z.print("Z:"); + arma::mat Z = arma::eye(4000,4000); +// 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; } \ No newline at end of file diff --git a/source/noise.c b/source/noise.c new file mode 100644 index 0000000..00f8d4f --- /dev/null +++ b/source/noise.c @@ -0,0 +1,80 @@ +// -------------------------------------------------------------- +// -------------------------------------------------------------- + +#include +#include +#include +#include + +#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; +} diff --git a/source/noise.h b/source/noise.h new file mode 100644 index 0000000..cc56e29 --- /dev/null +++ b/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_