first version

attic/eigen_demo.cpp

@@ -0,0 +1,95 @@
+/*******************************************************
+ A simple program that demonstrates the Eigen library.
+ The program defines a random symmetric matrix
+ and computes its eigendecomposition.
+ For further details read the Eigen Reference Manual
+********************************************************/
+
+#include <stdlib.h>
+#include <time.h>
+#include <string.h>
+
+// the following two are needed for printing
+#include <iostream>
+#include <iomanip>
+/**************************************
+/* The Eigen include files         */
+#include <Eigen/Core>
+#include <Eigen/Eigenvalues>
+#include <Eigen/QR>
+/***************************************/
+
+using namespace std;
+using namespace Eigen;
+
+int main(int argc, char **argv) {
+	int M = 3, N = 5;
+	MatrixXd X(M,N); // Define an M x N general matrix
+
+	// Fill X by random numbers between 0 and 9
+	// Note that indexing into matrices in NewMat is 1-based!
+	srand(time(NULL));
+	for (int i = 0; i < M; ++i) {
+		for (int j = 0; j < N; ++j) {
+			X(i,j) = rand() % 10;
+		}
+	}
+
+	MatrixXd C;
+	C = X * X.transpose(); // fill in C by X * X^t. 
+
+	cout << "The symmetrix matrix C" << endl;
+    cout << C << endl;
+
+
+	// compute eigendecomposition of C
+    SelfAdjointEigenSolver<MatrixXd> es(C);
+
+    MatrixXd D = es.eigenvalues().asDiagonal();
+    MatrixXd V = es.eigenvectors();
+
+	// Print the result
+	cout << "The eigenvalues matrix:" << endl;
+	cout << D << endl;
+	cout << "The eigenvectors matrix:" << endl;
+	cout << V << endl;
+
+	// Check that the first eigenvector indeed has the eigenvector property
+	VectorXd v1(3);
+	v1(0) = V(0,0);
+	v1(1) = V(1,0);
+	v1(2) = V(2,0);
+
+	VectorXd Cv1 = C * v1;
+	VectorXd lambda1_v1 = D(0) * v1;
+
+	cout << "The max-norm of the difference between C*v1 and lambda1*v1 is " << endl;
+	cout << Cv1.cwiseMax(lambda1_v1) << endl << endl;
+
+	// Build the inverse and check the result
+	MatrixXd Ci = C.inverse();
+	MatrixXd I  = Ci * C;
+
+	cout << "The inverse of C is" << endl;
+	cout << Ci << endl;
+	cout << "And the inverse times C is identity" << endl;
+	cout << I << endl;
+
+	// Example for multiple solves 
+	VectorXd r1(3), r2(3);
+	for (int i = 0; i < 3; ++i) {
+		r1(i) = rand() % 10;
+		r2(i) = rand() % 10;
+	}
+    ColPivHouseholderQR<MatrixXd> qr(C); // decomposes C
+	VectorXd s1 = qr.solve(r1);
+	VectorXd s2 = qr.solve(r2);
+
+	cout << "solution for right hand side r1" << endl;
+	cout << s1 << endl;
+	cout << "solution for right hand side r2" << endl;
+	cout << s2 << endl;
+
+	return 0;
+}
+

attic/eigen_map.cpp

@@ -0,0 +1,25 @@
+//
+// This example demonstrates the usage of Map in Eigen framework.
+// User-provided data is used to populate a Eigen matrix.
+//
+
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace Eigen;
+using namespace std;
+
+int main() {
+
+	// user data
+	double data[6] = { 1, 2, 3, 4, 5 , 6 };
+
+	// creating a 3x2 matrix from user data
+	MatrixXd mat = Map<MatrixXd>( data, 3, 2 );
+
+	// output it to see what is inside
+	cout << mat << endl;
+
+	return 0;
+}
+

attic/eigen_mult.cpp

@@ -0,0 +1,28 @@
+//
+// Example generates some random data in the Eigen matrix and outputs it.
+//
+
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace std;
+using namespace Eigen;
+
+int main() {
+
+	const int nCols = 5;
+	const int nRows = 4;
+
+	MatrixXd X( nRows, nCols );
+	X.setRandom();
+
+	MatrixXd Y( nCols, nRows );
+	Y.setRandom();
+
+	MatrixXd Z = X.matrix() * Y.matrix();
+
+	cout << Z << endl;
+
+	return 0;
+}
+

attic/eigen_random.cpp

@@ -0,0 +1,23 @@
+//
+// Example generates some random data in the Eigen matrix and outputs it.
+//
+
+#include <Eigen/Core>
+#include <iostream>
+
+using namespace std;
+using namespace Eigen;
+
+int main() {
+
+	const int nCols = 5;
+	const int nRows = 4;
+
+	MatrixXd X( nRows, nCols );
+	X.setRandom();
+
+	cout << X << endl;
+
+	return 0;
+}
+

attic/qsort_test.cpp

@@ -0,0 +1,33 @@
+#include <cmath>
+#include <iostream>
+
+using namespace std;
+
+// compare function
+int compare( const void *a, const void *b ) {
+
+	double eps = 0.00000001;
+	double diff = *(double *) a - *(double *) b;
+
+	if ( abs(diff) <= eps )
+		return 0;
+
+	if ( diff < 0 )
+		return -1;
+	else
+		return 1;
+}
+
+int main() {
+
+	double d[10] = { 0.1, 5.3, 4.2, 1.1, 0.9, 6.5, 0.1, 7.0, 5.9, 2.4 };
+
+	qsort( d, 10, sizeof(double), compare );
+
+	for ( int i = 0 ; i < 10 ; i++ )
+		cout << d[i] << " ";
+
+	cout << endl;
+
+	return 0;
+}

elml.hpp

@@ -0,0 +1,147 @@
+//
+// -----------------------------------------------
+// Provides Extreme Learning Machine algorithm.
+// -----------------------------------------------
+//
+// Refer to:
+// G.-B. Huang et al., “Extreme Learning Machine for Regression and Multiclass Classification,” IEEE TSMC, vol. 42, no. 2, pp. 513-529, 2012.
+//
+// Vladislavs D.
+//
+
+#include <iostream>
+#include <Eigen/Core>
+#include <Eigen/Cholesky>
+
+using namespace std;
+using namespace Eigen;
+
+#ifndef ELM_H
+#define ELM_H
+
+int compare( const void *a, const void *b );
+
+//template <typename Derived>
+MatrixXd buildTargetMatrix( double *Y, int nLabels );
+
+// entry function to train the ELM model
+// INPUT: X, Y, nhn, C
+// OUTPUT: inW, bias, outW
+template <typename Derived>
+int elmTrain( double *X, int dims, int nsmp, 
+       	      double *Y,
+	      const int nhn, const double C,
+	      MatrixBase<Derived> &inW, MatrixBase<Derived> &bias, MatrixBase<Derived> &outW ) {
+
+	// map the samples into the matrix object
+	MatrixXd mX = Map<MatrixXd>( X, dims, nsmp );
+
+	// build target matrix
+	MatrixXd mTargets = buildTargetMatrix( Y, nsmp );
+
+	// generate random input weight matrix - inW
+	inW = MatrixXd::Random( nhn, dims );
+
+	// generate random bias vectors
+	bias = MatrixXd::Random( nhn, 1 );
+
+	// compute the pre-H matrix
+	MatrixXd preH = inW * mX;
+
+	// compute hidden neuron output
+	MatrixXd H = (1 + (-preH.array()).exp()).cwiseInverse();
+
+	// build matrices to solve Ax = b
+	MatrixXd A = (MatrixXd::Identity(nhn,nhn)).array()*(1/C) + (H * H.transpose()).array();
+	MatrixXd b = H * mTargets.transpose();
+
+	// solve the output weights as a solution to a system of linear equations
+	outW = A.llt().solve( b );
+
+	return 0;
+
+}
+
+// entry function to predict class labels using the trained ELM model on test data
+// INPUT : X, inW, bias, outW
+// OUTPUT : scores
+template <typename Derived>
+int elmPredict( double *X, int dims, int nsmp,
+		MatrixBase<Derived> &mScores,
+		MatrixBase<Derived> &inW, MatrixBase<Derived> &bias, MatrixBase<Derived> &outW ) {
+
+	// map the sample into the Eigen's matrix object
+	MatrixXd mX = Map<MatrixXd>( X, dims, nsmp );
+
+	// build the pre-H matrix
+	MatrixXd preH = inW * mX + bias.replicate( 1, nsmp );
+
+	// apply the activation function
+	MatrixXd H = ( 1 + (-preH.array()).exp() ).cwiseInverse();
+
+	// compute output scores
+	mScores = (H.transpose() * outW).transpose();
+
+	return 0;
+}
+
+
+// --------------------------
+// Helper functions
+// --------------------------
+
+// compares two integer values
+//int compare( const void* a, const void *b ) {
+//	return ( *(int *) a - *(int *) b );
+//}
+
+int compare( const void *a, const void *b )
+{
+	const double *da = (const double *) a;
+	const double *db = (const double *) b;
+	return (*da > *db) - (*da < *db);
+}
+
+// builds 1-of-K target matrix from labels array
+//template <typename Derived>
+MatrixXd buildTargetMatrix( double *Y, int nLabels ) {
+
+	// make a temporary copy of the labels array
+	double *tmpY = new double[ nLabels ];
+	for ( int i = 0 ; i < nLabels ; i++ ) {
+		tmpY[i] = Y[i];
+	}
+
+	// sort the array of labels
+	qsort( tmpY, nLabels, sizeof(double), compare );
+
+
+	// count unique labels
+	int nunique = 1;
+	for ( int i = 0 ; i < nLabels - 1 ; i++ ) {
+		if ( tmpY[i] != tmpY[i+1] )
+			nunique++;
+	}
+
+	delete [] tmpY;
+
+	MatrixXd targets( nunique, nLabels );
+	targets.fill( 0 );
+
+
+	// fill in the ones
+	for ( int i = 0 ; i < nLabels ; i++ ) {
+		int idx = Y[i] - 1;
+		targets( idx, i ) = 1;
+	}
+
+	// normalize the targets matrix values (-1/1)
+	targets *= 2;
+	targets.array() -= 1;
+
+	return targets;
+
+}
+
+#endif
+