- 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

@@ -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)

RBM.jucer

@@ -20,9 +20,9 @@
     <LINUX_MAKE targetFolder="Builds/Linux">
       <CONFIGURATIONS>
         <CONFIGURATION name="Debug" libraryPath="/usr/X11R6/lib/" isDebug="1" optimisation="1"
-                       targetName="RBM"/>
+                       targetName="RBM" headerPath="/usr/local/include/eigen3"/>
         <CONFIGURATION name="Release" libraryPath="/usr/X11R6/lib/" isDebug="0" optimisation="3"
-                       targetName="RBM"/>
+                       targetName="RBM" headerPath="/usr/local/include/eigen3"/>
       </CONFIGURATIONS>
       <MODULEPATHS>
         <MODULEPATH id="juce_core" path="../../../../../../JUCE/modules"/>

Source/.cproject

@@ -22,9 +22,13 @@
 							<builder id="org.eclipse.cdt.build.core.settings.default.builder.1820544125" keepEnvironmentInBuildfile="false" managedBuildOn="false" name="Gnu Make Builder" superClass="org.eclipse.cdt.build.core.settings.default.builder"/>
 							<tool id="org.eclipse.cdt.build.core.settings.holder.libs.750397548" name="holder for library settings" superClass="org.eclipse.cdt.build.core.settings.holder.libs"/>
 							<tool id="org.eclipse.cdt.build.core.settings.holder.906730718" name="Assembly" superClass="org.eclipse.cdt.build.core.settings.holder">
+								<option id="org.eclipse.cdt.build.core.settings.holder.incpaths.852926707" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath"/>
 								<inputType id="org.eclipse.cdt.build.core.settings.holder.inType.1467109729" languageId="org.eclipse.cdt.core.assembly" languageName="Assembly" sourceContentType="org.eclipse.cdt.core.asmSource" superClass="org.eclipse.cdt.build.core.settings.holder.inType"/>
 							</tool>
 							<tool id="org.eclipse.cdt.build.core.settings.holder.1370151108" name="GNU C++" superClass="org.eclipse.cdt.build.core.settings.holder">
+								<option id="org.eclipse.cdt.build.core.settings.holder.incpaths.1839057859" name="Include Paths" superClass="org.eclipse.cdt.build.core.settings.holder.incpaths" valueType="includePath">
+									<listOptionValue builtIn="false" value="/usr/local/include/eigen3"/>
+								</option>
 								<inputType id="org.eclipse.cdt.build.core.settings.holder.inType.1834882069" languageId="org.eclipse.cdt.core.g++" languageName="GNU C++" sourceContentType="org.eclipse.cdt.core.cxxSource,org.eclipse.cdt.core.cxxHeader" superClass="org.eclipse.cdt.build.core.settings.holder.inType"/>
 							</tool>
 							<tool id="org.eclipse.cdt.build.core.settings.holder.280954803" name="GNU C" superClass="org.eclipse.cdt.build.core.settings.holder">

Source/.settings/language.settings.xml

@@ -5,8 +5,8 @@
 			<provider copy-of="extension" id="org.eclipse.cdt.ui.UserLanguageSettingsProvider"/>
 			<provider-reference id="org.eclipse.cdt.core.ReferencedProjectsLanguageSettingsProvider" ref="shared-provider"/>
 			<provider-reference id="org.eclipse.cdt.managedbuilder.core.MBSLanguageSettingsProvider" ref="shared-provider"/>
-			<provider copy-of="extension" id="org.eclipse.cdt.managedbuilder.core.GCCBuildCommandParser"/>
 			<provider-reference id="org.eclipse.cdt.managedbuilder.core.GCCBuiltinSpecsDetector" ref="shared-provider"/>
+			<provider copy-of="extension" id="org.eclipse.cdt.managedbuilder.core.GCCBuildCommandParser"/>
 		</extension>
 	</configuration>
 </project>

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];
+	m_data.resize(width*height);
-	memset(m_pData, 0, m_width*m_height*sizeof(double));
     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;
+    m_data.fill(0);
-    for (int i=0; i < m_width*m_height; i++)
+	setData(m_data);
-    {
-       	*(pData++) = 0;
-    }
-	setData(m_pData);
 }
 
-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<double>(std::max<double>(*pData, 0), 1);
+        	a = std::min<double>(std::max<double>((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();

Source/DrawComponent.h

@@ -30,6 +30,9 @@ public:
 	virtual ~DrawListener() {}
 	virtual void onDraw(DrawComponent &obj) = 0;
 };
+#include <Eigen/Dense>
+
+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);
+	void setData(const VectorXd& data);
-    const double* getData();
+	const VectorXd& getData();
 	void clear();
     //[/UserMethods]
 
@@ -79,7 +82,7 @@ private:
     float m_scaleX;
     float m_scaleY;
     ScopedPointer<Graphics>m_pG;
-    ScopedPointer<double>m_pData;
+    VectorXd m_data;
 	Image m_image;
 	double m_currData;
 	Colour m_currColor;

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;
-		double energy = 0;
+
+		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&)weights).hiddenBias()[index];
-		double sum = weights.getBiasVisible()[index];
-		const double *pStates = layer.getStates();
 
-		for (i=0; i <  layer.getNumUnits(); i++)
+		sum += layer.states().transpose() * ((Weights&)weights).weights().col(index);
-		{
-			sum += pStates[i] * weights.getWeights()[index][i];
-		}
 
 		return sum;
 	}

Source/Layer.hpp

@@ -11,19 +11,28 @@
 #ifndef LAYER_HPP
 #define LAYER_HPP
 #include <stdint.h>
+#include <iostream>
+#include <Eigen/Dense>
 #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_probs(numUnits)
-	, m_pStates(nullptr)
+	, 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;
+			return;
-			delete [] m_pStates;
 		}
 		m_numUnits = numUnits;
-		if (m_numUnits)
+		m_probs.resize(numUnits);
-		{
+		m_states.resize(numUnits);
-			m_pProbs = new double[m_numUnits];
-			m_pStates = new double[m_numUnits];
 
-			probsInit(0);
+		probsInit(0);
-			if (pStatesInit)
+		statesInit(0);
-			{
-				memcpy(m_pStates, pStatesInit, m_numUnits*sizeof(double));
-			}
-			else
-			{
-				statesInit(0);
-			}
-		}
 	}
 
-	Layer& operator= (const Layer &rhs)
+	void probsInit(const double &value)
 	{
-		memcpy(m_pProbs, rhs.m_pProbs, m_numUnits*sizeof(double));
+		m_probs.fill(value);
-		memcpy(m_pStates, rhs.m_pStates, m_numUnits*sizeof(double));
-
-		return *this;
 	}
 
-	Layer& operator+= (const Layer &rhs)
+	void statesInit(const double &value)
 	{
-		uint32_t i;
+		m_states.fill(value);
-		for (i=0; i < m_numUnits; i++)
-		{
-			m_pStates[i] += rhs.m_pStates[i];
-		}
-
-		return *this;
 	}
 
-	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:
+	inline double gaussProb(double x, double var) const
-	double *m_pProbs;
+	{
-	double *m_pStates;
+		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));
+	}
 
 
 };

Source/LayerArray.hpp

@@ -11,6 +11,9 @@
 #ifndef LAYERARRAY_HPP
 #define LAYERARRAY_HPP
 #include <stdint.h>
+#include <Eigen/Dense>
+
+using namespace Eigen;
 
 template <class T>
 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;
 		}
 

Source/MainComponent.cpp

@@ -21,6 +21,8 @@
 #ifdef WIN32
 #include <Windows.h>
 #endif
+#include <iostream>
+#include <Eigen/Dense>
 //[/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);
+        	pH = (const VectorXd*)&m_pRbm->toHidden(*pV);
-        	DrawHidden->setData(pH);
+        	DrawHidden->setData(*pH);
-        	pV = m_pRbm->toVisible(pH);
+        	pV = (const VectorXd*)&m_pRbm->toVisible(*pH);
-    		Draw2->setData(pV);
+    		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();
+		VectorXd w = m_weights.weights().col((int)sliderThatWasMoved->getValue());
-		if (!ppW)
+		VectorXd temp = w;
-			return;
-
-		double *pW = ppW[(int)sliderThatWasMoved->getValue()];
-		double *pTemp = new double [m_vNumX*m_vNumY];
 		double min = +1E12;
 		double max = -1E12;
 
 		for (int i=0; i < m_vNumX*m_vNumY; i++)
 		{
-    		pTemp[i] = pW[i];
+    		min = std::min<double>(min, (double)temp[i]);
-    		min = std::min<double>(min, pW[i]);
+    		max = std::max<double>(max, (double)temp[i]);
-    		max = std::max<double>(max, pW[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);
+		DrawWeights->setData(temp);
-		delete pTemp;
         //[/UserSliderCode_WeightsSlider]
     }
     else if (sliderThatWasMoved == numGibbsSlider)

Source/Rbm.hpp

@@ -12,6 +12,9 @@
 #include "HiddenLayer.hpp"
 #include "Weights.hpp"
 #include <cmath>
+#include <Eigen/Dense>
+
+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;
+		MatrixXd &w = (MatrixXd&)m_w.weights();
-		double dw;
-
-		double **ppW = m_w.getWeights();
-		const double *pH = h.getStates();
-		const double *pV = v.getStates();
 
 		// Update weights
-		for (i=0; i < m_w.getNumHidden(); i++)
+		w += mu*(v.states() * h.states().transpose());
-		{
+	}
-			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;
-		}
 
-		pBias = m_w.getBiasHidden();
+	void visibleBiasUpdate(VisibleLayer &v, double mu)
-		for (i=0; i < m_w.getNumHidden(); i++)
+	{
-		{
+		m_w.visibleBias().array() += mu*v.states().array();
-			dw = pH[i];
+	}
-			pBias[i] += mu*dw;
+
-		}
+	void hiddenBiasUpdate(HiddenLayer &h, double mu)
+	{
+		m_w.hiddenBias().array() += mu*h.states().array();
 	}
 
 	void train(LayerArray<VisibleLayer> &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++)
+			cout << h[j].probs() << endl;
-			{
-				p = h[j].getProbs()[i];
-				printf("%3.6f ", p);
-			}
-			printf("\n");
 		}
-		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++)
+			cout << h[j].states() << endl;
-			{
-				p = h[j].getStates()[i];
-				printf("%3.6f ", p);
-			}
-			printf("\n");
 		}
-		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++)
+			cout << vts.getAt(j).probs() << endl;
-			{
-				p = vts.getAt(j).getProbs()[i];
-				printf("%3.6f ", p);
-			}
-			printf("\n");
 		}
 
 		delete [] h;
 	}
 
-	const double* toHidden(const double *pVisible)
+	const VectorXd& toHidden(const VectorXd& visible)
 	{
-		double p;
+		VisibleLayer tv(m_w.getNumVisible(), (const VectorXd*)&visible);
-		uint32_t i;
-
-		VisibleLayer tv(m_w.getNumVisible(), pVisible);
 
 		m_th.probsUpdate(tv, m_w);
-		m_th.statesAssignfromProbs();
 
-#if 0
+		return m_th.probs();
-		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();
 	}
 
-	const double* toVisible(const double *pHidden)
+	const VectorXd& toVisible(const VectorXd& hidden)
 	{
-		double p;
+		HiddenLayer th(m_w.getNumHidden(), (const VectorXd*)&hidden);
-		uint32_t i;
-
-		HiddenLayer th(m_w.getNumHidden(), pHidden);
 
 		m_tv.probsUpdate(th, m_w);
-		m_tv.statesAssignfromProbs();
 
-#if 0
+		return m_tv.probs();
-		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();
 	}
 
 private:

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;
-		double energy = 0;
+
+		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&)weights).visibleBias()[index];
-		double sum = weights.getBiasVisible()[index];
-		const double *pStates = layer.getStates();
 
-		for (i=0; i <  layer.getNumUnits(); i++)
+		sum += layer.states().transpose() * ((Weights&)weights).weights().row(index).transpose();
-		{
-			sum += pStates[i] * weights.getWeights()[i][index];
-		}
 
 		return sum;
 	}

Source/Weights.hpp

@@ -11,48 +11,41 @@
 #ifndef WEIGHTS_HPP
 #define WEIGHTS_HPP
 #include <stdint.h>
+#include <iostream>
+#include <Eigen/Dense>
+
 #include "noise.h"
 
+using namespace std;
+using namespace Eigen;
+
 class Weights
 {
 public:
 	Weights(const char *pFilename)
-	: m_ppW(nullptr)
+	: m_numVisible(0)
-	, m_pBiasVisible(nullptr)
-	, m_pBiasHidden(nullptr)
-	, m_numVisible(0)
 	, m_numHidden(0)
 	{
 		load(pFilename);
 	}
 
 	Weights(uint32_t numVisible, uint32_t numHidden)
-	: m_ppW(nullptr)
+	: m_numVisible(numVisible)
-	, m_pBiasVisible(nullptr)
-	, m_pBiasHidden(nullptr)
-	, m_numVisible(numVisible)
 	, m_numHidden(numHidden)
 	{
 		Noise_Init(&m_noise, 0x32727155);
 		alloc(numVisible, numHidden);
-		shuffle(0);
 	}
 
 	Weights()
-	: m_ppW(nullptr)
+	: m_numVisible(0)
-	, m_pBiasVisible(nullptr)
-	, m_pBiasHidden(nullptr)
-	, m_numVisible(0)
 	, m_numHidden(0)
 	{
 		Noise_Init(&m_noise, 0x32727155);
 	}
 
 	Weights(const Weights &src)
-	: m_ppW(nullptr)
+	: m_numVisible(0)
-	, m_pBiasVisible(nullptr)
-	, m_pBiasHidden(nullptr)
-	, 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 ", m_w(i,j));
-				printf("%3.6f ", w);
 			}
 			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", m_bv(i));
-			printf("%3.6f\n", w);
 		}
 		printf("\n");
 
 		printf("bh = \n");
 		for (i=0; i < m_numHidden; i++)
 		{
-			w = m_pBiasHidden[i];
+			printf("%3.6f\n", m_bh(i));
-			printf("%3.6f\n", w);
 		}
 		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_numVisible == numVisible) && (m_numHidden == numHidden))
-
-		if (m_ppW)
 		{
-			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;
 	}
 };