interface4.cpp

#include <iostream>
#include <typeinfo>
#include <memory>
#include <algorithm>
#include <functional>
#include <random>
#include <Eigen/Eigen>

using namespace std;
using namespace Eigen;

/**
 * defines required typedefs and primitive type related things
 * which should be put under datatype and parameter code in Shogun
 */
typedef unsigned int index_t;
typedef int int32_t;
typedef float float32_t;
typedef double float64_t;

// these already exist in Shougn - see DataType.h
enum PTYPE
{
	INT32 = 0,
	FLOAT32,
	FLOAT64
};

/**
 * Dot product computers - all backend implementation are to be put here
 * can be put under shogun/mathematics/linalg/global/ or so...
 */
template <class T> class SGVector;

// dot product computers - base
template <class T, class Vector>
class VectorDotProduct
{
public:
	VectorDotProduct<T, Vector>() {}
	virtual T compute(Vector vector1, Vector vector2) const = 0;
	virtual ~VectorDotProduct<T, Vector>() {}
};

// eigen3 implementation
template <class T>
class Eigen3DotProduct : public VectorDotProduct<T, SGVector<T> >
{
public:
	Eigen3DotProduct<T>() {}
	virtual T compute(SGVector<T> vector1, SGVector<T> vector2) const
	{
		cout << "Eigen3DotProduct::compute() called" << endl;
		if (vector1.vlen != vector2.vlen)
			cerr << "dimension mismatch" << endl;
		Map<Matrix<T, Dynamic, 1> > vec1(vector1.vector.get(), vector1.vlen);
		Map<Matrix<T, Dynamic, 1> > vec2(vector2.vector.get(), vector2.vlen);
		return vec1.dot(vec2);
	}
	virtual ~Eigen3DotProduct<T>() {}
};

// naive implementation
template <class T>
class NaiveDotProduct : public VectorDotProduct<T, SGVector<T> >
{
public:
	NaiveDotProduct<T>() {}
	virtual T compute(SGVector<T> vector1, SGVector<T> vector2) const
	{
		cout << "NaiveDotProduct::compute() called" << endl;
		if (vector1.vlen != vector2.vlen)
			cerr << "dimension mismatch" << endl;
		T result = static_cast<T>(0);
		for (index_t i = 0; i < vector1.vlen; ++i)
			result += vector1[i] * vector2[i];
		return result;
	}
	virtual ~NaiveDotProduct<T>() {}
};

/**
 * global linear algebra backend setters
 */
enum ELinAlgBackend {Eigen3, Naive };

class LinearAlgebra
{
public:
	LinearAlgebra() {}
	// can also pass arguments for the datatypes for which we want to
	// initialize things - doesn't have to be for all datatypes
	void set_backend(ELinAlgBackend linalg_backend)
	{
		dot_computers = new void*[3]();
		switch (linalg_backend)
		{
			case Eigen3:
				dot_computers[INT32] = new Eigen3DotProduct<int32_t>();
				dot_computers[FLOAT32] = new Eigen3DotProduct<float32_t>();
				dot_computers[FLOAT64] = new Eigen3DotProduct<float64_t>();
				break;
			case Naive:
				dot_computers[INT32] = new NaiveDotProduct<int32_t>();
				dot_computers[FLOAT32] = new NaiveDotProduct<float32_t>();
				dot_computers[FLOAT64] = new NaiveDotProduct<float64_t>();
				break;
			default:
				break;
		}
	}
	template <class T, class Vector>
	VectorDotProduct<T, Vector>* get_dot_computer()
	{
	}
	~LinearAlgebra()
	{
		delete[] dot_computers;
	}
private:
	void **dot_computers;
};

#define GET_DOT_COMPUTER(T, Vector, PTYPE) \
template<> \
VectorDotProduct<T, Vector>* LinearAlgebra::get_dot_computer<T, Vector>() \
{ \
	return reinterpret_cast<VectorDotProduct<T, Vector>*>(dot_computers[PTYPE]); \
}
GET_DOT_COMPUTER(int32_t, SGVector<int32_t>, INT32)
GET_DOT_COMPUTER(float32_t, SGVector<float32_t>, FLOAT32)
GET_DOT_COMPUTER(float64_t, SGVector<float64_t>, FLOAT64)
#undef GET_DOT_COMPUTER

/** global linear algebra pointer - has to be inside init or so */
LinearAlgebra* sg_linalg;

/** the following should be modified inside Shogun to initialize */
void init_shogun()
{
	if (!sg_linalg)
		sg_linalg = new LinearAlgebra();
}

void exit_shogun()
{
	delete sg_linalg;
}

/**
 * mathematics/linalg modules
 */

// modified linear operator implementation
template <class Operand, class RetType>
class LinearOperator
{
public:
	virtual RetType apply(Operand operand) = 0;
};

// now vector dot operator - everything fits within apply then
template <class T, class Vector>
class VectorDotOperator : public LinearOperator<Vector, T>
{
public:
	VectorDotOperator(Vector vec) : vector(vec) {}
	virtual T apply(Vector other)
	{
		return sg_linalg->get_dot_computer<T, Vector>()->compute(vector, other);
	}
private:
	Vector vector;
};

/** shogun/lib/SGVector.h */
// could have used shogun's SGVector here instead - does nothing special now
template <class T>
class SGVector
{
public:
	SGVector<T>(index_t len) : vlen(len)
	{
		vector = shared_ptr<T>(new T[vlen](), [](T* ptr) { delete[] ptr; });
	}
	inline const T operator[](index_t idx) const { return vector.get()[idx]; }
	inline T operator[](index_t idx) { return vector.get()[idx]; }
	void set_const(T val)
	{
		for_each(&vector.get()[0], &vector.get()[vlen], [&val](T& elem) { elem = val; });
	}
	void display_vector(const char* name = "")
	{
		cout << name << ":";
		for_each(&vector.get()[0], &vector.get()[vlen], [](T& elem) { cout << elem << " "; });
		cout << endl;
	}
	virtual ~SGVector<T>() {}
	index_t vlen;
	shared_ptr<T> vector;
};

// have another implementation GPUVector or so...


/** now the main function */
int main()
{
	init_shogun();

	sg_linalg->set_backend(Eigen3);

	const index_t size = 10, iter = 10;
	SGVector<float64_t> v1(size), v2(size);
	default_random_engine generator;
	normal_distribution<float64_t> distrib(0.0, 1.0);

	VectorDotOperator<float64_t, SGVector<float64_t> > dotoperator(v1);

	for (index_t i = 0; i < iter; ++i)
	{
		for_each(&v1.vector.get()[0], &v1.vector.get()[size],
				[&distrib, &generator](float64_t& elem) { elem = distrib(generator); });
		for_each(&v2.vector.get()[0], &v2.vector.get()[size],
				[&distrib, &generator](float64_t& elem) { elem = distrib(generator); });
		v1.display_vector("v1");
		v2.display_vector("v2");
		cout << "v1.v2 : " << dotoperator.apply(v2) << endl;
	}

	exit_shogun();

	return 0;
}