From 19c69ac02adab231eaa14fa599f804fe4cabcd4a Mon Sep 17 00:00:00 2001 From: Jens Ahrensfeld Date: Sun, 12 Oct 2014 14:30:40 +0000 Subject: [PATCH] - use Matrix, linear algebra library Eigen 3.2.2 git-svn-id: http://moon:8086/svn/software/trunk/projects/RBM@23 b431acfa-c32f-4a4a-93f1-934dc6c82436 --- Builds/Linux/Makefile | 8 +- RBM.jucer | 4 +- Source/.cproject | 4 + Source/.settings/language.settings.xml | 2 +- Source/DrawComponent.cpp | 30 ++--- Source/DrawComponent.h | 9 +- Source/HiddenLayer.hpp | 22 ++-- Source/Layer.hpp | 126 ++++++++------------ Source/LayerArray.hpp | 23 ++-- Source/MainComponent.cpp | 44 ++++--- Source/Rbm.hpp | 159 +++++++++++-------------- Source/VisibleLayer.hpp | 22 ++-- Source/Weights.hpp | 152 ++++++++--------------- 13 files changed, 246 insertions(+), 359 deletions(-) diff --git a/Builds/Linux/Makefile b/Builds/Linux/Makefile index cfb00c0..0de0a10 100644 --- a/Builds/Linux/Makefile +++ b/Builds/Linux/Makefile @@ -18,12 +18,12 @@ ifeq ($(CONFIG),Debug) TARGET_ARCH := -march=native endif - CPPFLAGS := $(DEPFLAGS) -D "LINUX=1" -D "DEBUG=1" -D "_DEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules + CPPFLAGS := $(DEPFLAGS) -D "LINUX=1" -D "DEBUG=1" -D "_DEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules -I /usr/local/include/eigen3 CFLAGS += $(CPPFLAGS) $(TARGET_ARCH) -g -ggdb -O0 CXXFLAGS += $(CFLAGS) LDFLAGS += $(TARGET_ARCH) -L$(BINDIR) -L$(LIBDIR) -L/usr/X11R6/lib/ -lGL -lX11 -lXext -lXinerama -lasound -ldl -lfreetype -lpthread -lrt LDDEPS := - RESFLAGS := -D "LINUX=1" -D "DEBUG=1" -D "_DEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules + RESFLAGS := -D "LINUX=1" -D "DEBUG=1" -D "_DEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules -I /usr/local/include/eigen3 TARGET := RBM BLDCMD = $(CXX) -o $(OUTDIR)/$(TARGET) $(OBJECTS) $(LDFLAGS) $(RESOURCES) $(TARGET_ARCH) CLEANCMD = rm -rf $(OUTDIR)/$(TARGET) $(OBJDIR) @@ -39,12 +39,12 @@ ifeq ($(CONFIG),Release) TARGET_ARCH := -march=native endif - CPPFLAGS := $(DEPFLAGS) -D "LINUX=1" -D "NDEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules + CPPFLAGS := $(DEPFLAGS) -D "LINUX=1" -D "NDEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules -I /usr/local/include/eigen3 CFLAGS += $(CPPFLAGS) $(TARGET_ARCH) -O3 CXXFLAGS += $(CFLAGS) LDFLAGS += $(TARGET_ARCH) -L$(BINDIR) -L$(LIBDIR) -fvisibility=hidden -L/usr/X11R6/lib/ -lGL -lX11 -lXext -lXinerama -lasound -ldl -lfreetype -lpthread -lrt LDDEPS := - RESFLAGS := -D "LINUX=1" -D "NDEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules + RESFLAGS := -D "LINUX=1" -D "NDEBUG=1" -D "JUCER_LINUX_MAKE_7346DA2A=1" -D "JUCE_APP_VERSION=1.0.0" -D "JUCE_APP_VERSION_HEX=0x10000" -I /usr/include -I /usr/include/freetype2 -I ../../JuceLibraryCode -I ../../JuceLibraryCode/modules -I /usr/local/include/eigen3 TARGET := RBM BLDCMD = $(CXX) -o $(OUTDIR)/$(TARGET) $(OBJECTS) $(LDFLAGS) $(RESOURCES) $(TARGET_ARCH) CLEANCMD = rm -rf $(OUTDIR)/$(TARGET) $(OBJDIR) diff --git a/RBM.jucer b/RBM.jucer index 9cf50cb..70f0ece 100644 --- a/RBM.jucer +++ b/RBM.jucer @@ -20,9 +20,9 @@ + targetName="RBM" headerPath="/usr/local/include/eigen3"/> + targetName="RBM" headerPath="/usr/local/include/eigen3"/> diff --git a/Source/.cproject b/Source/.cproject index 71042e3..34721cb 100644 --- a/Source/.cproject +++ b/Source/.cproject @@ -22,9 +22,13 @@ + + diff --git a/Source/.settings/language.settings.xml b/Source/.settings/language.settings.xml index 980e504..89f540b 100644 --- a/Source/.settings/language.settings.xml +++ b/Source/.settings/language.settings.xml @@ -5,8 +5,8 @@ - + diff --git a/Source/DrawComponent.cpp b/Source/DrawComponent.cpp index 567871c..2682a41 100644 --- a/Source/DrawComponent.cpp +++ b/Source/DrawComponent.cpp @@ -48,11 +48,11 @@ DrawComponent::DrawComponent (int width, int height) //[Constructor] You can add your own custom stuff here.. // m_pG = new Graphics(m_image); - m_pData = new double[width*height]; - memset(m_pData, 0, m_width*m_height*sizeof(double)); + m_data.resize(width*height); Noise_Init(&noise, 0x3231); m_pG->setImageResamplingQuality(Graphics::lowResamplingQuality); + clear(); //[/Constructor] } @@ -65,7 +65,6 @@ DrawComponent::~DrawComponent() //[Destructor]. You can add your own custom destruction code here.. m_pG = nullptr; - m_pData = nullptr; //[/Destructor] } @@ -195,7 +194,7 @@ void DrawComponent::drawAt(int x, int y, bool setColor) { if (setColor) { - if (m_pData[index] > 0.0) + if (m_data[index] > 0.0) { m_currData = 0.0; m_currColor = (Colours::black); @@ -206,7 +205,7 @@ void DrawComponent::drawAt(int x, int y, bool setColor) m_currColor = (Colours::white); } } - m_pData[index] = m_currData; + m_data[index] = m_currData; m_pG->setColour (m_currColor); } } @@ -219,35 +218,26 @@ void DrawComponent::drawAt(int x, int y, bool setColor) void DrawComponent::clear() { - double *pData = m_pData; - for (int i=0; i < m_width*m_height; i++) - { - *(pData++) = 0; - } - setData(m_pData); + m_data.fill(0); + setData(m_data); } -const double* DrawComponent::getData () +const VectorXd& DrawComponent::getData () { - return m_pData; + return m_data; } -void DrawComponent::setData (const double *pData) +void DrawComponent::setData (const VectorXd& data) { double a; - for (int i=0; i < m_width*m_height; i++) - { - m_pData[i] = pData[i]; - } for (int i=0; i < m_height; i++) { for (int j=0; j < m_width; j++) { - a = std::min(std::max(*pData, 0), 1); + a = std::min(std::max((double)data[i*m_width + j], 0), 1); m_pG->setColour(Colour(Colours::white).greyLevel(a)); m_pG->fillRect(m_scaleX*j, m_scaleY*i, m_scaleX, m_scaleY); - pData++; } } repaint(); diff --git a/Source/DrawComponent.h b/Source/DrawComponent.h index 1c0d6a4..b23ac98 100644 --- a/Source/DrawComponent.h +++ b/Source/DrawComponent.h @@ -30,6 +30,9 @@ public: virtual ~DrawListener() {} virtual void onDraw(DrawComponent &obj) = 0; }; +#include + +using namespace Eigen; //[/Headers] @@ -53,8 +56,8 @@ public: //[UserMethods] -- You can add your own custom methods in this section. void setListener(DrawListener *pListener); void drawAt(int x, int y, bool setColor); - void setData(const double *pData); - const double* getData(); + void setData(const VectorXd& data); + const VectorXd& getData(); void clear(); //[/UserMethods] @@ -79,7 +82,7 @@ private: float m_scaleX; float m_scaleY; ScopedPointerm_pG; - ScopedPointerm_pData; + VectorXd m_data; Image m_image; double m_currData; Colour m_currColor; diff --git a/Source/HiddenLayer.hpp b/Source/HiddenLayer.hpp index 29bc182..649ee9a 100644 --- a/Source/HiddenLayer.hpp +++ b/Source/HiddenLayer.hpp @@ -14,7 +14,7 @@ class HiddenLayer : public Layer { public: - HiddenLayer(uint32_t numUnits = 0, const double *pStatesInit = nullptr) + HiddenLayer(uint32_t numUnits = 0, const VectorXd *pStatesInit = nullptr) : Layer(numUnits, pStatesInit) { } @@ -25,27 +25,19 @@ public: double getEnergy(const Weights &weights) { - uint32_t i; - double energy = 0; + double energy; + + energy = -((Weights&)weights).hiddenBias().transpose() * states(); - for (i=0; i < getNumUnits(); i++) - { - energy -= weights.getBiasHidden()[i] * getStates()[i]; - } return energy; } private: - double accum(const Layer &layer, const Weights &weights, uint32_t index) const + double accum(Layer &layer, Weights &weights, uint32_t index) { - uint32_t i; - double sum = weights.getBiasVisible()[index]; - const double *pStates = layer.getStates(); + double sum = ((Weights&)weights).hiddenBias()[index]; - for (i=0; i < layer.getNumUnits(); i++) - { - sum += pStates[i] * weights.getWeights()[index][i]; - } + sum += layer.states().transpose() * ((Weights&)weights).weights().col(index); return sum; } diff --git a/Source/Layer.hpp b/Source/Layer.hpp index a952138..7dee41a 100644 --- a/Source/Layer.hpp +++ b/Source/Layer.hpp @@ -11,19 +11,28 @@ #ifndef LAYER_HPP #define LAYER_HPP #include +#include +#include #include "noise.h" #include "Weights.hpp" +using namespace Eigen; + class Layer { public: - Layer(uint32_t numUnits = 0, const double *pStatesInit = nullptr) + Layer(uint32_t numUnits = 0, const VectorXd *pStatesInit = nullptr) : m_numUnits(numUnits) - , m_pProbs(nullptr) - , m_pStates(nullptr) + , m_probs(numUnits) + , m_states(numUnits) { - setNumUnits(numUnits, pStatesInit); + setNumUnits(numUnits); Noise_Init(&m_noise, 0x12345677); + + if (pStatesInit && (pStatesInit->size() == numUnits)) + { + m_states = *pStatesInit; + } } virtual ~Layer() @@ -32,95 +41,55 @@ public: Noise_Free(&m_noise); } - void setNumUnits(uint32_t numUnits, const double *pStatesInit = nullptr) + void setNumUnits(uint32_t numUnits) { - if (m_numUnits) + if (m_numUnits == numUnits) { - delete [] m_pProbs; - delete [] m_pStates; + return; } m_numUnits = numUnits; - if (m_numUnits) - { - m_pProbs = new double[m_numUnits]; - m_pStates = new double[m_numUnits]; + m_probs.resize(numUnits); + m_states.resize(numUnits); - probsInit(0); - if (pStatesInit) - { - memcpy(m_pStates, pStatesInit, m_numUnits*sizeof(double)); - } - else - { - statesInit(0); - } - } + probsInit(0); + statesInit(0); } - Layer& operator= (const Layer &rhs) + void probsInit(const double &value) { - memcpy(m_pProbs, rhs.m_pProbs, m_numUnits*sizeof(double)); - memcpy(m_pStates, rhs.m_pStates, m_numUnits*sizeof(double)); - - return *this; + m_probs.fill(value); } - Layer& operator+= (const Layer &rhs) + void statesInit(const double &value) { - uint32_t i; - for (i=0; i < m_numUnits; i++) - { - m_pStates[i] += rhs.m_pStates[i]; - } - - return *this; + m_states.fill(value); } - void probsInit(double value) const + void probsUpdate(Layer &layer, Weights &weights, double lambda = 1.0, double variance = 1.0) + { + probsUpdateLogistic(layer, weights, lambda, variance); + } + + void probsUpdateLogistic(Layer &layer, Weights &weights, double lambda, double variance) { uint32_t i; for (i=0; i < m_numUnits; i++) { - m_pProbs[i] = value; + m_probs(i) = logSigmoid(lambda/variance*accum(layer, weights, i)); } } - void statesInit(double value) const + void probsUpdateGaussian(Layer &layer, Weights &weights, double lambda, double variance) { uint32_t i; for (i=0; i < m_numUnits; i++) { - m_pStates[i] = value; + m_probs(i) = gaussProb(lambda*accum(layer, weights, i), variance); } } - void probsUpdate(const Layer &layer, const Weights &weights) const - { - uint32_t i; - - for (i=0; i < m_numUnits; i++) - { - m_pProbs[i] = logSigmoid(accum(layer, weights, i)); - } - } - - void statesScale(double kscale) const - { - uint32_t i; - - for (i=0; i < m_numUnits; i++) - { - m_pStates[i] *= kscale; - } - } - - void statesAssignfromProbs() - { - memcpy(m_pStates, m_pProbs, m_numUnits*sizeof(double)); - } - void statesUpdateStochastic() { uint32_t i; @@ -129,18 +98,18 @@ public: for (i=0; i < m_numUnits; i++) { sample = Noise_Uniform(&m_noise, 0.5); - m_pStates[i] = (double)(sample <= m_pProbs[i]); + m_states(i) = (double)(sample <= m_probs(i)); } } - const double *getProbs() const + VectorXd& probs() { - return m_pProbs; + return m_probs; } - const double *getStates() const + VectorXd& states() { - return m_pStates; + return m_states; } uint32_t getNumUnits() const @@ -151,19 +120,28 @@ public: virtual double getEnergy(const Weights &weights) = 0; private: - uint32_t m_numUnits; noise_gen_t m_noise; +protected: + uint32_t m_numUnits; + VectorXd m_probs; + VectorXd m_states; + + virtual double accum(Layer &layer, Weights &weights, uint32_t index) = 0; inline double logSigmoid(double x) const { return 1./(1 + exp(-x)); } -protected: - double *m_pProbs; - double *m_pStates; + inline double gaussProb(double x, double var) const + { + double k = 1.0/sqrt(var*2*3.14159265359); - virtual double accum(const Layer &layer, const Weights &weights, uint32_t index) const = 0; + double mu = 0; + double x2 = (x-mu)*(x-mu); + + return k*exp(-x2/(2*var)); + } }; diff --git a/Source/LayerArray.hpp b/Source/LayerArray.hpp index 9e899a5..4fb4832 100644 --- a/Source/LayerArray.hpp +++ b/Source/LayerArray.hpp @@ -11,6 +11,9 @@ #ifndef LAYERARRAY_HPP #define LAYERARRAY_HPP #include +#include + +using namespace Eigen; template class LayerArray; @@ -31,7 +34,7 @@ class LayerArray class Entry : public T { public: - Entry(uint32_t numUnits, const double *pInit) + Entry(uint32_t numUnits, const VectorXd *pInit) : T(numUnits, pInit) , pPrev(nullptr) , pNext(nullptr) @@ -75,13 +78,15 @@ public: clear(); } - void add(const double *pData, uint32_t size) + T* add(const VectorXd *pData, uint32_t size) { + Entry *pNew; Entry *pL = m_pRoot; + pNew = new Entry(size, pData); if (!m_pRoot) { - m_pRoot = new Entry(size, pData); + m_pRoot = pNew; } else { @@ -90,12 +95,14 @@ public: { pL = pL->pNext; } - pL->pNext = new Entry(size, pData); + pL->pNext = pNew; pL->pNext->pPrev = pL; } rebuildIndex(); if (m_pListener) m_pListener->onChanged(*this); + + return (T*)pNew; } void clear() @@ -172,7 +179,7 @@ public: for (j=0; j < numUnits; j++) { - fprintf(pFile, "%3.6f\n", m_ppIndex[i]->getStates()[j]); + fprintf(pFile, "%3.6f\n", m_ppIndex[i]->states()(j)); } } @@ -199,13 +206,15 @@ public: uint32_t numUnits; fscanf(pFile, "%d\n", &numUnits); pData = new double[numUnits]; + T* data = add(nullptr, numUnits); + for (j=0; j < numUnits; j++) { fscanf(pFile, "%f", &v); - pData[j] = v; + data->states()(j) = v; } - add(pData, numUnits); + delete [] pData; } diff --git a/Source/MainComponent.cpp b/Source/MainComponent.cpp index 59d926d..a56e680 100644 --- a/Source/MainComponent.cpp +++ b/Source/MainComponent.cpp @@ -21,6 +21,8 @@ #ifdef WIN32 #include #endif +#include +#include //[/Headers] #include "MainComponent.h" @@ -50,6 +52,8 @@ void mylog(const char* format, ...) //[/MiscUserDefs] +using namespace Eigen; +using namespace std; //============================================================================== MainComponent::MainComponent () : m_layers(this), @@ -220,12 +224,12 @@ MainComponent::MainComponent () projectNameLabel->setText(String("TestPrj"), dontSendNotification ); numEpochslabel->setText(String(100), dontSendNotification ); - learningRateLabel->setText(String(0.2), dontSendNotification ); + learningRateLabel->setText(String(0.1), dontSendNotification ); rbmUseExpectationsToggleButton->setToggleState(false, sendNotification); rbmDoRaoBlackwellToggleButton->setToggleState(false, sendNotification); rbmDoRobbinsMonroToggleButton->setToggleState(false, sendNotification); rbmReduceVarianceToggleButton->setToggleState(false, sendNotification); - //[/Constructor] + //[/Constructor] } MainComponent::~MainComponent() @@ -334,7 +338,7 @@ void MainComponent::buttonClicked (Button* buttonThatWasClicked) { //[UserButtonCode_addButton] -- add your button handler code here.. Draw2->setData(Draw->getData()); - m_layers.add(Draw->getData(), m_vNumX*m_vNumY); + m_layers.add(&Draw->getData(), m_vNumX*m_vNumY); patterSlider->setRange (0, m_layers.getSize()-1, 1); //[/UserButtonCode_addButton] } @@ -411,15 +415,15 @@ void MainComponent::buttonClicked (Button* buttonThatWasClicked) { //[UserButtonCode_reconstructEquButton] -- add your button handler code here.. uint32_t i; - const double *pV, *pH; + const VectorXd *pV, *pH; - pV = Draw->getData(); + pV = (const VectorXd*)&Draw->getData(); for (i=0; i < 100; i++) { - pH = m_pRbm->toHidden(pV); - DrawHidden->setData(pH); - pV = m_pRbm->toVisible(pH); - Draw2->setData(pV); + pH = (const VectorXd*)&m_pRbm->toHidden(*pV); + DrawHidden->setData(*pH); + pV = (const VectorXd*)&m_pRbm->toVisible(*pH); + Draw2->setData(*pV); } //[/UserButtonCode_reconstructEquButton] } @@ -463,39 +467,33 @@ void MainComponent::sliderValueChanged (Slider* sliderThatWasMoved) if (m_layers.getSize() > 0) { VisibleLayer &p = (VisibleLayer&)m_layers.getAt((int)sliderThatWasMoved->getValue()); - Draw2->setData(p.getStates()); + Draw2->setData(p.states()); } //[/UserSliderCode_patterSlider] } else if (sliderThatWasMoved == WeightsSlider) { //[UserSliderCode_WeightsSlider] -- add your slider handling code here.. - double **ppW = m_weights.getWeights(); - if (!ppW) - return; - - double *pW = ppW[(int)sliderThatWasMoved->getValue()]; - double *pTemp = new double [m_vNumX*m_vNumY]; + VectorXd w = m_weights.weights().col((int)sliderThatWasMoved->getValue()); + VectorXd temp = w; double min = +1E12; double max = -1E12; for (int i=0; i < m_vNumX*m_vNumY; i++) { - pTemp[i] = pW[i]; - min = std::min(min, pW[i]); - max = std::max(max, pW[i]); + min = std::min(min, (double)temp[i]); + max = std::max(max, (double)temp[i]); } for (int i=0; i < m_vNumX*m_vNumY; i++) { - pTemp[i] -= min; + temp[i] -= min; } for (int i=0; i < m_vNumX*m_vNumY; i++) { - pTemp[i] /= (max-min); + temp[i] /= (max-min); } - DrawWeights->setData(pTemp); - delete pTemp; + DrawWeights->setData(temp); //[/UserSliderCode_WeightsSlider] } else if (sliderThatWasMoved == numGibbsSlider) diff --git a/Source/Rbm.hpp b/Source/Rbm.hpp index ba3ca23..96ce223 100644 --- a/Source/Rbm.hpp +++ b/Source/Rbm.hpp @@ -12,6 +12,9 @@ #include "HiddenLayer.hpp" #include "Weights.hpp" #include +#include + +using namespace Eigen; void mylog(const char* format, ...); #define printf mylog @@ -47,35 +50,20 @@ public: void weightsUpdate(VisibleLayer &v, HiddenLayer &h, double mu) { - uint32_t i, j; - double dw; - - double **ppW = m_w.getWeights(); - const double *pH = h.getStates(); - const double *pV = v.getStates(); + MatrixXd &w = (MatrixXd&)m_w.weights(); // Update weights - for (i=0; i < m_w.getNumHidden(); i++) - { - for (j=0; j < m_w.getNumVisible(); j++) - { - dw = pV[j] * pH[i]; - ppW[i][j] += mu*dw; - } - } - double *pBias = m_w.getBiasVisible(); - for (i=0; i < m_w.getNumVisible(); i++) - { - dw = pV[i]; - pBias[i] += mu*dw; - } + w += mu*(v.states() * h.states().transpose()); + } - pBias = m_w.getBiasHidden(); - for (i=0; i < m_w.getNumHidden(); i++) - { - dw = pH[i]; - pBias[i] += mu*dw; - } + void visibleBiasUpdate(VisibleLayer &v, double mu) + { + m_w.visibleBias().array() += mu*v.states().array(); + } + + void hiddenBiasUpdate(HiddenLayer &h, double mu) + { + m_w.hiddenBias().array() += mu*h.states().array(); } void train(LayerArray &vt, uint32_t numEpochs, double mu, uint32_t numGibbs = 1, bool useExpectations = false, bool doRaoBlackwell = false, bool useProbsForHiddenReconstruction = false, bool doRobbinsMonro = false) @@ -93,6 +81,10 @@ public: double dProgress = 1.0/numEpochs; m_progress = 0; + const double lambda = 1.0; + const double variance = 1.0; + const double penalty = 0.0; + if (useExpectations) { mu /= vt.getSize(); @@ -103,7 +95,7 @@ public: for (i=0; i < vt.getSize(); i++) { // Create hidden layer base on training data - ht[i].probsUpdate(vt[i], w); + ht[i].probsUpdate(vt[i], w, lambda, variance); } } @@ -112,7 +104,7 @@ public: for (i=0; i < vt.getSize(); i++) { t = (uint32_t)(0.5 + (vt.getSize()-1)*Noise_Uniform(&m_noise, 0.5)); - h.probsUpdate(vt[t], w); + h.probsUpdate(vt[t], w, lambda, variance); // Create hidden layer base on training data if (doRobbinsMonro) @@ -127,23 +119,25 @@ public: // Update weights (positive phase) if (doRaoBlackwell) { - pH->statesAssignfromProbs(); + pH->states() = pH->probs(); } else { pH->statesUpdateStochastic(); } weightsUpdate(vt[t], h, +mu); + visibleBiasUpdate(vt[t], +mu); + hiddenBiasUpdate(h, +mu); for (gibbs=0; gibbs < numGibbs; gibbs++) { pH->statesUpdateStochastic(); // Create visible reconstruction (a fantasy...) - v.probsUpdate(*pH, w); + v.probsUpdate(*pH, w, lambda, variance); if (useProbsForHiddenReconstruction) { - v.statesAssignfromProbs(); + v.states() = v.probs(); } else { @@ -151,18 +145,20 @@ public: } // Create hidden reconstruction - pH->probsUpdate(v, w); + pH->probsUpdate(v, w, lambda, variance); } // Update weights (negative phase) if (doRaoBlackwell) { - pH->statesAssignfromProbs(); + pH->states() = pH->probs(); } else { pH->statesUpdateStochastic(); } weightsUpdate(v, *pH, -mu); + visibleBiasUpdate(v, -mu); + hiddenBiasUpdate(*pH, -mu); if (!useExpectations) { @@ -170,11 +166,28 @@ public: } } +#if 0 + { + HiddenLayer th(m_w.getNumHidden()); + for (i=0; i < vt.getSize(); i++) + { + ht[i].probsUpdate(vt[t], w, lambda, variance); + ht[i].statesUpdateStochastic(); + th += ht[i]; + } + th *= 1.0/vt.getSize(); + th += -0.02; + hiddenBiasUpdate(th, -mu); + } +#endif + if (useExpectations) { w = m_w; } + getEnergy(v, *pH); + m_progress += dProgress; if (m_pListener) { @@ -199,9 +212,12 @@ public: { for (j=0; j < v.getNumUnits(); j++) { - energy -= v.getStates()[j] * h.getStates()[i] * m_w.getWeights()[i][j]; +// energy -= v.getStates()[j] * h.getStates()[i] * m_w.getWeights()[i][j]; } } + + // ToDo: make this correct +// energy -= (v.states().transpose() * h.states()); // * m_w.weights(); return energy; } @@ -223,28 +239,18 @@ public: } printf("pi(t) = (pi^, v>)\n"); - for (i=0; i < m_w.getNumHidden(); i++) + for (j=0; j < vts.getSize(); j++) { - for (j=0; j < vts.getSize(); j++) - { - p = h[j].getProbs()[i]; - printf("%3.6f ", p); - } - printf("\n"); + cout << h[j].probs() << endl; } - printf("\n"); + cout << endl; printf("si(t) = (si^, v>)\n"); - for (i=0; i < m_w.getNumHidden(); i++) + for (j=0; j < vts.getSize(); j++) { - for (j=0; j < vts.getSize(); j++) - { - p = h[j].getStates()[i]; - printf("%3.6f ", p); - } - printf("\n"); + cout << h[j].states() << endl; } - printf("\n"); + cout << endl; printf("p(v) = (t^, v>)\n"); for (i=0; i < vts.getSize(); i++) @@ -257,11 +263,11 @@ public: for (j=0; j < vts.getSize(); j++) { p = exp(-getEnergy(vts.getAt(j), h[i]))/z; - printf("%3.6f ", p); + cout << p << endl; } - printf("\n"); + cout << endl; } - printf("\n"); + cout << endl; // Reconstruct for (i=0; i < vts.getSize(); i++) @@ -270,61 +276,30 @@ public: } printf("A fantasy... (v^, t>)\n"); - for (i=0; i < m_w.getNumVisible(); i++) + for (j=0; j < vts.getSize(); j++) { - for (j=0; j < vts.getSize(); j++) - { - p = vts.getAt(j).getProbs()[i]; - printf("%3.6f ", p); - } - printf("\n"); + cout << vts.getAt(j).probs() << endl; } delete [] h; } - const double* toHidden(const double *pVisible) + const VectorXd& toHidden(const VectorXd& visible) { - double p; - uint32_t i; - - VisibleLayer tv(m_w.getNumVisible(), pVisible); + VisibleLayer tv(m_w.getNumVisible(), (const VectorXd*)&visible); m_th.probsUpdate(tv, m_w); - m_th.statesAssignfromProbs(); -#if 0 - printf("pi(t) = (pi^, v>)\n"); - for (i=0; i < m_w.getNumHidden(); i++) - { - p = m_th.getProbs()[i]; - printf("%3.6f\n", p); - } - printf("\n"); -#endif - return m_th.getStates(); + return m_th.probs(); } - const double* toVisible(const double *pHidden) + const VectorXd& toVisible(const VectorXd& hidden) { - double p; - uint32_t i; - - HiddenLayer th(m_w.getNumHidden(), pHidden); + HiddenLayer th(m_w.getNumHidden(), (const VectorXd*)&hidden); m_tv.probsUpdate(th, m_w); - m_tv.statesAssignfromProbs(); -#if 0 - printf("pi(t) = (pi^, v>)\n"); - for (i=0; i < m_w.getNumVisible(); i++) - { - p = m_tv.getProbs()[i]; - printf("%3.6f\n", p); - } - printf("\n"); -#endif - return m_tv.getStates(); + return m_tv.probs(); } private: diff --git a/Source/VisibleLayer.hpp b/Source/VisibleLayer.hpp index 502ae42..99ee593 100644 --- a/Source/VisibleLayer.hpp +++ b/Source/VisibleLayer.hpp @@ -14,7 +14,7 @@ class VisibleLayer : public Layer { public: - VisibleLayer(uint32_t numUnits = 0, const double *pStatesInit = nullptr) + VisibleLayer(uint32_t numUnits = 0, const VectorXd *pStatesInit = nullptr) : Layer(numUnits, pStatesInit) { } @@ -25,27 +25,19 @@ public: double getEnergy(const Weights &weights) { - uint32_t i; - double energy = 0; + double energy; + + energy = -((Weights&)weights).visibleBias().transpose() * states(); - for (i=0; i < getNumUnits(); i++) - { - energy -= weights.getBiasVisible()[i] * getStates()[i]; - } return energy; } private: - double accum(const Layer &layer, const Weights &weights, uint32_t index) const + double accum(Layer &layer, Weights &weights, uint32_t index) { - uint32_t i; - double sum = weights.getBiasVisible()[index]; - const double *pStates = layer.getStates(); + double sum = ((Weights&)weights).visibleBias()[index]; - for (i=0; i < layer.getNumUnits(); i++) - { - sum += pStates[i] * weights.getWeights()[i][index]; - } + sum += layer.states().transpose() * ((Weights&)weights).weights().row(index).transpose(); return sum; } diff --git a/Source/Weights.hpp b/Source/Weights.hpp index c638f6e..b869c0d 100644 --- a/Source/Weights.hpp +++ b/Source/Weights.hpp @@ -11,48 +11,41 @@ #ifndef WEIGHTS_HPP #define WEIGHTS_HPP #include +#include +#include + #include "noise.h" +using namespace std; +using namespace Eigen; + class Weights { public: Weights(const char *pFilename) - : m_ppW(nullptr) - , m_pBiasVisible(nullptr) - , m_pBiasHidden(nullptr) - , m_numVisible(0) + : m_numVisible(0) , m_numHidden(0) { load(pFilename); } Weights(uint32_t numVisible, uint32_t numHidden) - : m_ppW(nullptr) - , m_pBiasVisible(nullptr) - , m_pBiasHidden(nullptr) - , m_numVisible(numVisible) + : m_numVisible(numVisible) , m_numHidden(numHidden) { Noise_Init(&m_noise, 0x32727155); alloc(numVisible, numHidden); - shuffle(0); } Weights() - : m_ppW(nullptr) - , m_pBiasVisible(nullptr) - , m_pBiasHidden(nullptr) - , m_numVisible(0) + : m_numVisible(0) , m_numHidden(0) { Noise_Init(&m_noise, 0x32727155); } Weights(const Weights &src) - : m_ppW(nullptr) - , m_pBiasVisible(nullptr) - , m_pBiasHidden(nullptr) - , m_numVisible(0) + : m_numVisible(0) , m_numHidden(0) { Noise_Init(&m_noise, 0x32727155); @@ -69,7 +62,6 @@ public: void setUnits(uint32_t numVisible, uint32_t numHidden) { alloc(numVisible, numHidden); - shuffle(0); } void shuffle(double stdDev) @@ -77,68 +69,52 @@ public: uint32_t i, j; double kdev = stdDev*sqrt(12.0); - for (j=0; j < m_numVisible; j++) + for (i=0; i < m_numVisible; i++) { - m_pBiasVisible[j] = kdev*Noise_Uniform(&m_noise, 0.5); + m_bv(i) = kdev*Noise_Uniform(&m_noise, 0.5); } - for (i=0; i < m_numHidden; i++) + for (j=0; j < m_numHidden; j++) { - m_pBiasHidden[i] = kdev*Noise_Uniform(&m_noise, 0.5); + m_bh(j) = kdev*Noise_Uniform(&m_noise, 0.5); } - for (i=0; i < m_numHidden; i++) + for (i=0; i < m_numVisible; i++) { - for (j=0; j < m_numVisible; j++) + for (j=0; j < m_numHidden; j++) { - m_ppW[i][j] = kdev*Noise_Uniform(&m_noise, 0.5); + m_w(i,j) = kdev*Noise_Uniform(&m_noise, 0.5); } } } Weights& operator= (const Weights &rhs) { - uint32_t i, j; + m_bv = rhs.m_bv; + m_bh = rhs.m_bh; + m_w = rhs.m_w; - for (j=0; j < m_numVisible; j++) - { - m_pBiasVisible[j] = rhs.m_pBiasVisible[j]; - } - - for (i=0; i < m_numHidden; i++) - { - m_pBiasHidden[i] = rhs.m_pBiasHidden[i]; - } - - for (i=0; i < m_numHidden; i++) - { - for (j=0; j < m_numVisible; j++) - { - m_ppW[i][j] = rhs.m_ppW[i][j]; - } - } return *this; } - double **getWeights() const + MatrixXd& weights() { - return m_ppW; + return m_w; } - double *getBiasVisible() const + VectorXd& visibleBias() { - return m_pBiasVisible; + return m_bv; } - double *getBiasHidden() const + VectorXd& hiddenBias() { - return m_pBiasHidden; + return m_bh; } void print() { uint32_t i, j; - double w; printf("\n"); printf("w(v,h) = (v^, h>)\n"); @@ -146,8 +122,7 @@ public: { for (j=0; j < m_numHidden; j++) { - w = m_ppW[j][i]; - printf("%3.6f ", w); + printf("%3.6f ", m_w(i,j)); } printf("\n"); } @@ -156,18 +131,18 @@ public: printf("bv = \n"); for (i=0; i < m_numVisible; i++) { - w = m_pBiasVisible[i]; - printf("%3.6f\n", w); + printf("%3.6f\n", m_bv(i)); } printf("\n"); printf("bh = \n"); for (i=0; i < m_numHidden; i++) { - w = m_pBiasHidden[i]; - printf("%3.6f\n", w); + printf("%3.6f\n", m_bh(i)); } printf("\n"); + + cout << m_w << endl; } uint32_t getNumVisible() @@ -196,17 +171,17 @@ public: for (i=0; i < m_numVisible; i++) { - fprintf(pFile, "%3.6f\n", m_pBiasVisible[i]); + fprintf(pFile, "%3.6f\n", m_bv(i)); } for (i=0; i < m_numHidden; i++) { - fprintf(pFile, "%3.6f\n", m_pBiasHidden[i]); + fprintf(pFile, "%3.6f\n", m_bh(i)); } for (i=0; i < m_numVisible; i++) { for (j=0; j < m_numHidden; j++) { - fprintf(pFile, "%3.6f ", m_ppW[j][i]); + fprintf(pFile, "%3.6f ", m_w(i,j)); } fprintf(pFile, "\n"); } @@ -234,12 +209,12 @@ public: for (i=0; i < numVisible; i++) { fscanf(pFile, "%f", &v); - m_pBiasVisible[i] = v; + m_bv(i) = v; } for (i=0; i < numHidden; i++) { fscanf(pFile, "%f", &v); - m_pBiasHidden[i] = v; + m_bh(i) = v; } for (i=0; i < numVisible; i++) { @@ -247,7 +222,7 @@ public: { fscanf(pFile, "%f", &v); - m_ppW[j][i] = v; + m_w(i, j) = v; } } @@ -255,60 +230,31 @@ public: } private: - double **m_ppW; - double *m_pBiasVisible; - double *m_pBiasHidden; uint32_t m_numVisible; uint32_t m_numHidden; noise_gen_t m_noise; + MatrixXd m_w; + VectorXd m_bv; + VectorXd m_bh; void alloc(uint32_t numVisible, uint32_t numHidden) { - uint32_t i; - - if (m_ppW) + if ((m_numVisible == numVisible) && (m_numHidden == numHidden)) { - if ((numVisible == m_numVisible) && (numHidden == m_numHidden)) - { return; - } - free(); - alloc(numVisible, numHidden); - } - else - { - m_numVisible = numVisible; - m_numHidden = numHidden; - m_ppW = new double*[m_numHidden]; - for (i=0; i < m_numHidden; i++) - { - m_ppW[i] = new double[m_numVisible]; - } - m_pBiasVisible = new double[m_numVisible]; - m_pBiasHidden = new double[m_numHidden]; } + m_numVisible = numVisible; + m_numHidden = numHidden; + + m_w.resize(numVisible, numHidden); + m_bv.resize(numVisible); + m_bh.resize(numHidden); + + shuffle(0); } void free() { - uint32_t i; - - - if (m_ppW) - { - for (i=0; i < m_numHidden; i++) - { - delete [] m_ppW[i]; - } - delete [] m_ppW; - m_ppW = nullptr; - } - - delete [] m_pBiasVisible; - m_pBiasVisible = nullptr; - - delete [] m_pBiasHidden; - m_pBiasHidden = nullptr; } };