dot_factory_benchmark.cpp

#include <iostream>
#include <typeinfo>
#include <memory>
#include <algorithm>
#include <functional>
#include <random>
#include <Eigen/Eigen>
#include <shogun/base/init.h>
#include <shogun/lib/SGVector.h>
#include <shogun/lib/Time.h>
#include <shogun/io/SGIO.h>

using namespace std;
using namespace Eigen;
using namespace shogun;

/**
 * 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;
*/
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, vector1.vlen);
		Map<Matrix<T, Dynamic, 1> > vec2(vector2.vector, 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()
	{
		for (index_t i=0; i<3; ++i)
			delete[] dot_computers[i];
		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 my_init_shogun()
{
	if (!sg_linalg)
		sg_linalg = new LinearAlgebra();
}

void my_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()
{
	my_init_shogun();
	init_shogun_with_defaults();

	sg_linalg->set_backend(Eigen3);

	const index_t size = 1000000, iter = 1000;
	SGVector<float32_t> a(size), b(size);
	/*
	default_random_engine generator;
	normal_distribution<float32_t> distrib(0.0, 1.0);

	VectorDotOperator<float64_t, SGVector<float64_t> > dotoperator(v1);
*/
	CTime *time = new CTime();
	double mean = 0, var = 0;

	for (index_t i = 1; i <= iter; ++i)
	{
		SGVector<float>::random_vector(a.vector, a.vlen, 0, 1);
		SGVector<float>::random_vector(b.vector, b.vlen, 0, 1);
		/*
		for_each(&v1.vector[0], &v1.vector[size],
				[&distrib, &generator](float32_t& elem) { elem = distrib(generator); });
		for_each(&v2.vector[0], &v2.vector[size],
				[&distrib, &generator](float32_t& elem) { elem = distrib(generator); });
		v1.display_vector("v1");
		v2.display_vector("v2");
		cout << "v1.v2 : " << dotoperator.apply(v2) << endl;
		*/
		time->start();
		sg_linalg->get_dot_computer<float32_t, SGVector<float32_t> >()->compute(a, b);
		double elapsed = time->cur_time_diff();
		double delta = elapsed - mean;
		mean += delta/i;
		var += delta * (elapsed - mean);
	}
	var /= iter;
	SG_SPRINT("mean %.15f, var %.15lf\n", mean, var);

	my_exit_shogun();
	exit_shogun();

	return 0;
}

dot_macro_benchmark.cpp

#include <iostream>
#include <shogun/lib/SGVector.h>
#include <Eigen/Eigen>
#include <shogun/lib/Time.h>
#include <shogun/base/init.h>
#include <shogun/io/SGIO.h>

using namespace std;
using namespace shogun;
using namespace Eigen;

enum linalg { Naive, Eigen3 };

#ifndef LINALG_BACKEND
#define LINALG_BACKEND(Backend) \
namespace Use_##Backend \
{ \
	const linalg backend = Backend; \
}
#endif // LINALG_BACKEND

LINALG_BACKEND(Naive)
LINALG_BACKEND(Eigen3)

#undef LINALG_BACKEND

namespace shogun
{
class impl
{
	template <class Scalar, class Vector, enum linalg>
	struct dot_compute
	{
		static Scalar compute(Vector a, Vector b)
		{
			// cout << "Generic dot_compute" << endl;
			return static_cast<Scalar>(0);
		}
	};

	template <class Scalar>
	struct dot_compute<Scalar, SGVector<Scalar>, Eigen3>
	{
		static Scalar compute(SGVector<Scalar> a, SGVector<Scalar> b)
		{
			// cout << "Eigen3 dot_compute" << endl;
			typedef Matrix<Scalar, Dynamic, 1> VectorXt;
			Map<VectorXt> vec_a(a.vector, a.vlen);
			Map<VectorXt> vec_b(b.vector, b.vlen);
			return vec_a.dot(vec_b);
		}
	};
	friend class LinearAlgegra;
};
}

#ifndef SET_LINEAR_ALGEBRA_BACKEND
#define SET_LINEAR_ALGEBRA_BACKEND(BACKEND) \
namespace LinearAlgebra \
{ \
	using Use_##BACKEND::backend; \
}
#endif

#ifdef USE_EIGEN3
SET_LINEAR_ALGEBRA_BACKEND(Eigen3)
#else
SET_LINEAR_ALGEBRA_BACKEND(Naive)
#endif

namespace LinearAlgebra
{
	template <class Scalar, class Vector>
	static Scalar dot(Vector a, Vector b)
	{
		return impl::dot_compute<Scalar, Vector, backend>::compute(a, b);
	}
}

int main()
{
	init_shogun_with_defaults();
	index_t size = 1000000;
	SGVector<float> a(size), b(size);
	index_t iter = 1000;
	CTime* time = new CTime();
	double mean = 0, var = 0;
	for (index_t i = 1; i <= iter; ++i)
	{
		SGVector<float>::random_vector(a.vector, a.vlen, 0, 1);
		SGVector<float>::random_vector(b.vector, b.vlen, 0, 1);
		time->start();
		LinearAlgebra::dot<float, SGVector<float> >(a, b);
		double elapsed = time->cur_time_diff();
		double delta = elapsed - mean;
		mean += delta/i;
		var += delta * (elapsed - mean);
	}
	var /= iter;
	SG_SPRINT("mean %.15f, var %.15lf\n", mean, var);
	SG_UNREF(time);
	exit_shogun();
	return 0;
}