Merge pull request #1173 from lambday/feature/log_determinant

class CLinearOperator and DenseMatrixOperator added

src/shogun/mathematics/logdet/DenseMatrixOperator.cpp

@@ -0,0 +1,110 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * Written (W) 2013 Soumyajit De
+ */
+
+#include <shogun/lib/common.h>
+#include <shogun/mathematics/logdet/DenseMatrixOperator.h>
+#include <shogun/lib/SGVector.h>
+#include <shogun/lib/SGMatrix.h>
+
+#ifdef HAVE_EIGEN3
+#include <shogun/mathematics/eigen3.h>
+
+using namespace Eigen;
+#endif // HAVE_EIGEN3
+
+namespace shogun
+{
+
+template<class T>
+CDenseMatrixOperator<T>::CDenseMatrixOperator()
+	: CLinearOperator<T>()
+	{
+		SG_SGCDEBUG("%s created (%p)\n", this->get_name(), this);
+	}
+
+template<class T>
+CDenseMatrixOperator<T>::CDenseMatrixOperator(SGMatrix<T> op)
+	: CLinearOperator<T>(op.num_cols), m_operator(op)
+	{
+		SG_SGCDEBUG("%s created (%p)\n", this->get_name(), this);
+	}
+
+template<class T>
+CDenseMatrixOperator<T>::~CDenseMatrixOperator()
+	{
+		SG_SGCDEBUG("%s destroyed (%p)\n", this->get_name(), this);
+	}
+
+template<class T>
+SGMatrix<T> CDenseMatrixOperator<T>::get_matrix_operator() const
+	{
+		return m_operator;
+	}
+
+#ifdef HAVE_EIGEN3
+template<class T>
+SGVector<T> CDenseMatrixOperator<T>::apply(SGVector<T> b) const
+	{
+		REQUIRE(m_operator.matrix, "Operator not initialized!\n");
+		REQUIRE(this->get_dim()==b.vlen,
+			"Number of rows of vector must be equal to the "
+			"number of cols of the operator!\n");
+
+		typedef Matrix<T, Dynamic, 1> VectorXt;
+		typedef Matrix<T, Dynamic, Dynamic> MatrixXt;
+
+		Map<VectorXt> _b(b.vector, b.vlen);
+		Map<MatrixXt> _op(m_operator.matrix, m_operator.num_rows,
+			m_operator.num_cols);
+
+		SGVector<T> result(b.vlen);
+		Map<VectorXt> _result(result.vector, result.vlen);
+		_result=_op*_b;
+
+		return result;
+	}
+
+#define UNDEFINED(type) \
+template<> \
+SGVector<type> CDenseMatrixOperator<type>::apply(SGVector<type> b) const \
+	{	\
+		SG_SERROR("Not supported for %s\n", #type);\
+		return b; \
+	}
+
+UNDEFINED(bool)
+UNDEFINED(char)
+UNDEFINED(int8_t)
+UNDEFINED(uint8_t)
+#undef UNDEFINED
+
+#else
+template<class T>
+SGVector<T> CDenseMatrixOperator<T>::apply(SGVector<T> b) const
+	{
+		SG_SWARNING("Eigen3 required!\n");
+		return b;
+	}
+#endif // HAVE_EIGEN3
+
+template class CDenseMatrixOperator<bool>;
+template class CDenseMatrixOperator<char>;
+template class CDenseMatrixOperator<int8_t>;
+template class CDenseMatrixOperator<uint8_t>;
+template class CDenseMatrixOperator<int16_t>;
+template class CDenseMatrixOperator<uint16_t>;
+template class CDenseMatrixOperator<int32_t>;
+template class CDenseMatrixOperator<uint32_t>;
+template class CDenseMatrixOperator<int64_t>;
+template class CDenseMatrixOperator<uint64_t>;
+template class CDenseMatrixOperator<float32_t>;
+template class CDenseMatrixOperator<float64_t>;
+template class CDenseMatrixOperator<floatmax_t>;
+template class CDenseMatrixOperator<complex64_t>;
+}

src/shogun/mathematics/logdet/DenseMatrixOperator.h

@@ -0,0 +1,56 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * Written (W) 2013 Soumyajit De
+ */
+
+#ifndef DENSE_MATRIX_OPERATOR_H_
+#define DENSE_MATRIX_OPERATOR_H_
+
+#include <shogun/lib/config.h>
+#include <shogun/mathematics/logdet/LinearOperator.h>
+
+namespace shogun
+{
+template<class T> class SGVector;
+template<class T> class SGMatrix;
+
+/** @brief Class that represents a dense-matrix linear operator */
+template<class T> class CDenseMatrixOperator : public CLinearOperator<T>
+{
+public:
+	/** default constructor, no args */
+	CDenseMatrixOperator();
+
+	/** default constructor, one args */
+	CDenseMatrixOperator(SGMatrix<T> op);
+
+	/** destructor */
+	~CDenseMatrixOperator();
+
+	/** method that applies the dense-matrix linear operator to a vector
+	 * @param b the vector to which the linear operator applies
+	 * @return the result vector
+	 */
+	virtual SGVector<T> apply(SGVector<T> b) const;
+
+	/** @return the dense matrix operator */
+	SGMatrix<T> get_matrix_operator() const;
+
+	/** @return object name */
+	virtual const char* get_name() const
+	{
+		return "CDenseMatrixOperator";
+	}
+
+private:
+	/** the dense matrix operator */
+	const SGMatrix<T> m_operator;
+};
+
+}
+
+#endif // LINEAR_OPERATOR_H_

src/shogun/mathematics/logdet/LinearOperator.h

@@ -0,0 +1,90 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * Written (W) 2013 Soumyajit De
+ */
+#ifndef LINEAR_OPERATOR_H_
+#define LINEAR_OPERATOR_H_
+
+#include <shogun/lib/config.h>
+#include <shogun/base/SGObject.h>
+
+namespace shogun
+{
+template<class T> class SGVector;
+
+/** @brief Abstract template base class that represents a linear operator,
+ *  e.g. a matrix
+ */
+template<class T> class CLinearOperator : public CSGObject
+{
+public:
+	/** default constructor, no args */
+	CLinearOperator()
+	: CSGObject(), dim(0)
+	{
+		SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
+	}
+
+	/** default constructor, one arg */
+	CLinearOperator(index_t _dim)
+	: CSGObject(), dim(_dim)
+	{
+		SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
+	}
+
+	/** copy constructor */
+	CLinearOperator(const CLinearOperator<T>& orig)
+	: CSGObject(orig), dim(orig.get_dim())
+	{
+		SG_GCDEBUG("%s created (%p)\n", this->get_name(), this)
+	}
+
+	/** destructor */
+	virtual ~CLinearOperator()
+	{
+		SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this)
+	}
+
+	/** getter for the dimension */
+	const index_t get_dim() const
+	{
+		return dim;
+	}
+
+	/** abstract method that applies the linear operator to a vector
+	 * @param b the vector to which the linear operator applies
+	 * @return the result vector
+	 */
+	virtual SGVector<T> apply(SGVector<T> b) const = 0;
+
+	/** @return object name */
+	virtual const char* get_name() const
+	{
+		return "CLinearOperator";
+	}
+protected:
+	/** the dimension of vector on which the linear operator can apply */
+	const index_t dim;
+};
+
+template class CLinearOperator<bool>;
+template class CLinearOperator<char>;
+template class CLinearOperator<int8_t>;
+template class CLinearOperator<uint8_t>;
+template class CLinearOperator<int16_t>;
+template class CLinearOperator<uint16_t>;
+template class CLinearOperator<int32_t>;
+template class CLinearOperator<uint32_t>;
+template class CLinearOperator<int64_t>;
+template class CLinearOperator<uint64_t>;
+template class CLinearOperator<float32_t>;
+template class CLinearOperator<float64_t>;
+template class CLinearOperator<floatmax_t>;
+template class CLinearOperator<complex64_t>;
+}
+
+#endif // LINEAR_OPERATOR_H_

tests/unit/mathematics/logdet/DenseMatrixOperator_unittest.cc

@@ -0,0 +1,82 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * Written (W) 2013 Soumyajit De
+ */
+
+#include <shogun/lib/config.h>
+#include <shogun/lib/SGMatrix.h>
+#include <shogun/lib/SGVector.h>
+#include <shogun/mathematics/eigen3.h>
+#include <shogun/mathematics/logdet/DenseMatrixOperator.h>
+#include <gtest/gtest.h>
+
+using namespace shogun;
+
+#ifdef HAVE_EIGEN3
+using namespace Eigen;
+
+TEST(DenseMatrixOperator, apply)
+{
+	const index_t size=5;
+	SGVector<float64_t> b1(size);
+	b1.set_const(0.25);
+
+	// float64_t, fixed matrix
+	SGMatrix<float64_t> m1(size, size);
+	m1.set_const(0.5);
+
+	CDenseMatrixOperator<float64_t> op11(m1);
+	SGVector<float64_t> r1=op11.apply(b1);
+
+	for (index_t i=0; i<r1.vlen; ++i)
+	{
+		EXPECT_NEAR(r1[i], 0.625, 1E-16);
+	}
+
+	// float64_t, identity matrix
+	typedef Matrix<float64_t, size, size> MatrixSd;
+	Map<MatrixSd> mapped1(m1.matrix, m1.num_rows, m1.num_cols);
+	mapped1=MatrixSd::Identity();
+	CDenseMatrixOperator<float64_t> op12(m1);
+	r1=op12.apply(b1);
+
+	for (index_t i=0; i<r1.vlen; ++i)
+	{
+		EXPECT_NEAR(r1[i], 0.25, 1E-16);
+	}
+
+	SGVector<complex64_t> b2(size);
+	b2.set_const(complex64_t(0.25, 1.0));
+
+	// complex64_t, fixed matrix	
+	SGMatrix<complex64_t> m2(size, size);
+	m2.set_const(complex64_t(0.5, 0.25));
+
+	CDenseMatrixOperator<complex64_t> op21(m2);
+	SGVector<complex64_t> r2=op21.apply(b2);
+
+	for (index_t i=0; i<r2.vlen; ++i)
+	{
+		EXPECT_NEAR(r2[i].real(), -0.625, 1E-16);
+		EXPECT_NEAR(r2[i].imag(), 2.8125, 1E-16);
+	}
+
+	// complex64_t, identity matrix
+	typedef Matrix<complex64_t, size, size> MatrixScd;
+	Map<MatrixScd> mapped2(m2.matrix, m2.num_rows, m2.num_cols);
+	mapped2=MatrixScd::Identity();
+	CDenseMatrixOperator<complex64_t> op22(m2);
+	r2=op22.apply(b2);
+
+	for (index_t i=0; i<r2.vlen; ++i)
+	{
+		EXPECT_NEAR(r2[i].real(), 0.25, 1E-16);
+		EXPECT_NEAR(r2[i].imag(), 1.0, 1E-16);
+	}
+}
+
+#endif // HAVE_EIGEN3