objects.h

// -*- c++ -*-

// Copyright (C) 2005,2009 Tom Drummond (twd20@cam.ac.uk),
// Ed Rosten (er258@cam.ac.uk), Gerhard Reitmayr (gr281@cam.ac.uk)
//
// This file is part of the TooN Library.  This library 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 2, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License along
// with this library; see the file COPYING.  If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
// USA.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the GNU General Public License.

namespace TooN {

namespace Internal{
	// dummy structs that are used in 0-ary operators
	struct Zero;
	struct SizedZero;
	struct RCZero;
	template<class P> struct Identity;
	template<class P> struct SizedIdentity;

	template<int S, class P, class B, class Ps> class ScalarsVector;	
	template<int R, int C, class P, class B, class Ps> class ScalarsMatrix;	
	template<int R, int C, class P, class B, class Ps> class AddIdentity;	
	template<class P> class Scalars;	
	template<class P> class SizedScalars;	
	template<class P> class RCScalars;
	
	///@internal
	///@brief This class represents 1 and only in all its forms.
	///@ingroup gInternal
	struct One{

		One(){} ///<This constructor does nothing. This allows const One to be declared with no initializer.
		///Generic cast to anything
		template<class C> operator C() const
		{
			return 1;
		}
	};
	template<class Rhs> Rhs operator*(One, const Rhs& v){return v;}   ///<Multiplies One by something.
	template<class Lhs> Lhs operator*(const Lhs& v, One){return v;}   ///<Multiplies something by One
	template<class Rhs> Rhs operator+(One, const Rhs& v){return 1+v;} ///<Adds something to One
	template<class Lhs> Lhs operator+(const Lhs& v, One){return v+1;} ///<Adds One to something
	template<class Rhs> Rhs operator-(One, const Rhs& v){return 1-v;} ///<Subtracts something from One
	template<class Lhs> Lhs operator-(const Lhs& v, One){return v-1;} ///<Subtracts One from something.

	///Returns negative One.
	inline int operator-(const One&)
	{
		return -1;
	}
	
	///@internal
	///@brief For an instance \e i of type C, what is the type of \e -i?
	///Usually the answer is that is it the same type.
	///@ingroup gInternal
	template<class C> struct NegType
	{
		typedef C Type; ///<The type of -C
	};

	/**@internal
	   @brief The type of -One
	   @ingroup gInternal
	*/
	template<> struct NegType<One>
	{
		///One really repersents 1. Therefore -One is not the same type
		///as One.
		typedef int Type;
	};

	///@internal
	///@brief Does One behave as a field with respect to Rhs?
	///@ingroup gInternal
	template<class Rhs> struct Field<One, Rhs>
	{
		///One can be converted in to anything, so the resulting type is
		///a field if the other type is a field.
		static const int is = IsField<Rhs>::value;
	};

	///@internal
	///@brief Does One behave as a field with respect to Lhs?
	///@ingroup gInternal
	template<class Lhs> struct Field<Lhs, One>
	{
		///One can be converted in to anything, so the resulting type is
		///a field if the other type is a field.
		static const int is = IsField<Lhs>::value;
	};

}

////////////////////
// Zero
////////////////////



template<> struct Operator<Internal::SizedZero>;
template<> struct Operator<Internal::RCZero>;


///@internal
///@brief Object which behaves like a block of zeros. See TooN::Zeros.
///@ingroup gInternal
template<> struct Operator<Internal::Zero> {
	///@name Operator members
	///@{
	template<int Size, class Precision, class Base>
	void eval(Vector<Size, Precision, Base>& v) const {
		for(int i=0; i < v.size(); i++) {
			v[i]= 0;
		}
	}

	template<int R, int C, class P, class B>
	void eval(Matrix<R,C,P,B>& m) const {
		for(int r=0; r<m.num_rows(); r++){
			for(int c=0; c<m.num_cols(); c++){
				m(r,c)=0;
			}
		}
	}
	///@}

	template<int R, int C, class P, class B>
	bool notequal(Matrix<R,C,P,B>& m) const {
		for(int r=0; r<m.num_rows(); r++)
			for(int c=0; c<m.num_cols(); c++)
				if(m[r][c] != 0)
					return 1;
		
		return 0;
	}


	template<int S, class P, class B>
	bool notequal(Vector<S,P,B>& v) const {
		for(int i=0; i<v.size(); i++)
			if(v[i] != 0)
				return 1;
		return 0;
	}

	///Generate a sized Zero object for constructing dynamic vectors.
	Operator<Internal::SizedZero> operator()(int s);
	///Generate a sized Zero object for constructing dynamic matrices.
	Operator<Internal::RCZero> operator()(int r, int c);
	
};

///@internal
///@brief Variant of the Zeros object which holds two sizes for constructing dynamic matrices.
///@ingroup gInternal
template<> struct Operator<Internal::RCZero> : public Operator<Internal::Zero> {

	///@name Operator members determining the size.
	///@{
	Operator(int r, int c) : my_rows(r), my_cols(c) {}

	const int my_rows;
	const int my_cols;
	
	int num_rows() const {return my_rows;}
	int num_cols() const {return my_cols;}
	///@}
};

///@internal
///@brief Variant of the Zeros object which holds a size for constructing dynamic vectors.
///@ingroup gInternal
template<> struct Operator<Internal::SizedZero> : public Operator<Internal::Zero> {

	///@name Operator members determining the size for vectors and square matrices.
	///@{
	Operator(int s)	: my_size(s) {}
		
	const int my_size;
	
	int size() const {return my_size;}
	int num_rows() const {return my_size;}
	int num_cols() const {return my_size;}
	///@}
};

inline Operator<Internal::SizedZero> Operator<Internal::Zero>::operator()(int s){
	return Operator<Internal::SizedZero>(s);
}

inline Operator<Internal::RCZero> Operator<Internal::Zero>::operator()(int r, int c){
	return Operator<Internal::RCZero>(r,c);
}


//////////////
// Identity
//////////////

///@internal
///@brief Operator to construct a new matrix with idendity added 
///@ingroup gInternal
template<int R, int C, class P, class B, class Precision> struct Operator<Internal::AddIdentity<R,C,P,B,Precision> >
{
	const Precision s;   ///<Scale of the identity matrix
	const Matrix<R,C,P,B>& m; ///<matrix to which the identity should be added
	bool invert_m; ///<Whether the identity should be added to + or - m
	
	///@name Construction
	///@{
	Operator(Precision s_, const Matrix<R,C,P,B>& m_, bool b)
		:s(s_),m(m_),invert_m(b){}
	///@}

	///@name Operator members
	///@{
	template<int R1, int C1, class P1, class B1>
	void eval(Matrix<R1,C1,P1,B1>& mm) const{
		for(int r=0; r < m.num_rows(); r++)
			for(int c=0; c < m.num_cols(); c++)
				if(invert_m)
					mm[r][c] = -m[r][c];
				else
					mm[r][c] = m[r][c];

		for(int i=0; i < m.num_rows(); i++)
				mm[i][i] += (P)s;
	}
	///@}

	///@name Sized operator members
	///@{
	int num_rows() const
	{
		return m.num_rows();
	}
	int num_cols() const
	{
		return m.num_cols();
	}
	///@}
};

///@internal
///@brief Object which behaves like an Identity matrix. See TooN::Identity.
///@ingroup gInternal
template<class Pr> struct Operator<Internal::Identity<Pr> > {
	
	///@name Scalable operators members
	///@{

	typedef Pr Precision;
	template<class Pout, class Pmult> Operator<Internal::Identity<Pout> > scale_me(const Pmult& m) const
	{
		return Operator<Internal::Identity<Pout> >(val*m);
	}
	///}
	
	///<Scale of the identity matrix.
	const Precision val;

	///@name Construction
	///@{
	Operator(const Precision& v)
		:val(v)
	{}

	Operator()
	{}
	///}

	///@name Operator members
	///@{
	template<int R, int C, class P, class B>
	void eval(Matrix<R,C,P,B>& m) const {
		SizeMismatch<R, C>::test(m.num_rows(), m.num_cols());
		
		for(int r=0; r<m.num_rows(); r++){
			for(int c=0; c<m.num_cols(); c++){
				m(r,c)=0;
			}
		}
		
		for(int r=0; r < m.num_rows(); r++) {
			m(r,r) = (P)val;
		}
	}
	
	template<int Rows, int Cols, typename P, typename B>
	void plusequals(Matrix<Rows, Cols, P, B>& m) const
	{
		SizeMismatch<Rows, Cols>::test(m.num_rows(), m.num_cols());
		for(int i=0; i < m.num_rows(); i++)
			m[i][i] += (P)val;
	}

	template <int Rows, int Cols, typename P1, typename B1> 
	Operator<Internal::AddIdentity<Rows,Cols,P1,B1,Precision> > add(const Matrix<Rows,Cols, P1, B1>& m) const
	{
		SizeMismatch<Rows, Cols>::test(m.num_rows(), m.num_cols());
		return Operator<Internal::AddIdentity<Rows,Cols,P1,B1,Precision> >(val, m, 0);
	}

	template <int Rows, int Cols, typename P1, typename B1> 
	Operator<Internal::AddIdentity<Rows,Cols,P1,B1,Precision> > rsubtract(const Matrix<Rows,Cols, P1, B1>& m) const
	{
		SizeMismatch<Rows, Cols>::test(m.num_rows(), m.num_cols());
		return Operator<Internal::AddIdentity<Rows,Cols,P1,B1,Precision> >(-val, m, 0);
	}

	template <int Rows, int Cols, typename P1, typename B1> 
	Operator<Internal::AddIdentity<Rows,Cols,P1,B1,Precision> > lsubtract(const Matrix<Rows,Cols, P1, B1>& m) const
	{
		SizeMismatch<Rows, Cols>::test(m.num_rows(), m.num_cols());
		return Operator<Internal::AddIdentity<Rows,Cols,P1,B1,Precision> >(val, m, 1);
	}
	///@}
	
	///@name Sizeable operator members
	///@{
	Operator<Internal::SizedIdentity<Precision> > operator()(int s){
		return Operator<Internal::SizedIdentity<Precision> >(s, val);
	}
	///@}
};
	
///@internal
///@brief A variant of Identity which holds a size, allowing dynamic matrices to be constructed
///@ingroup gInternal
template<class Precision> struct Operator<Internal::SizedIdentity<Precision> > 
	: public  Operator<Internal::Identity<Precision> > {

	using Operator<Internal::Identity<Precision> >::val;
	
	///@name Constructors
	///@{
	Operator(int s, const Precision& v)
		:Operator<Internal::Identity<Precision> > (v), my_size(s)
	{}
	///@}

	///@name Sized operator members
	///@{
	const int my_size;
	int num_rows() const {return my_size;}
	int num_cols() const {return my_size;}
	///@}

	///@name Scalable operator members
	///@{
	template<class Pout, class Pmult> Operator<Internal::SizedIdentity<Pout> > scale_me(const Pmult& m) const
	{
		return Operator<Internal::SizedIdentity<Pout> >(my_size, val*m);
	}
	///@}
};
////////////////////////////////////////////////////////////////////////////////
//
// Addition of scalars to vectors and matrices
//

	
///@internal
///@brief Operator to construct a new vector a a vector with a scalar added to every element
///@ingroup gInternal
template<int S, class P, class B, class Precision> struct Operator<Internal::ScalarsVector<S,P,B,Precision> >
{
	const Precision s;      ///<Scalar to add
	const Vector<S,P,B>& v; ///<Vector to be added to.
	const bool invert_v;    ///<Whether to use + or - \c v

	///@name Constructors
	///@{
	Operator(Precision s_, const Vector<S,P,B>& v_, bool inv)
		:s(s_),v(v_),invert_v(inv){}
	///@}

	///@name Operator members
	///@{
	template<int S1, class P1, class B1>
	void eval(Vector<S1,P1,B1>& vv) const{
		for(int i=0; i < v.size(); i++)
			if(invert_v)
				vv[i] = s - v[i];
			else
				vv[i] = s + v[i];
	}
	///@}

	///@name Sized operator members
	///@{
	int size() const
	{
		return v.size();
	}
	///@}
};

///@internal
///@brief Operator to construct a new matrix a a matrix with a scalar added to every element
///@ingroup gInternal
template<int R, int C, class P, class B, class Precision> struct Operator<Internal::ScalarsMatrix<R,C,P,B,Precision> >
{
	const Precision s;        ///<Scalar to add
	const Matrix<R,C,P,B>& m; ///<Vector to be added to.
	const bool invert_m;      ///<Whether to use + or - \c m
	///@name Operator members
	///@{
	Operator(Precision s_, const Matrix<R,C,P,B>& m_, bool inv)
		:s(s_),m(m_),invert_m(inv){}
	template<int R1, int C1, class P1, class B1>
	void eval(Matrix<R1,C1,P1,B1>& mm) const{
		for(int r=0; r < m.num_rows(); r++)
			for(int c=0; c < m.num_cols(); c++)
				if(invert_m)
					mm[r][c] = s - m[r][c];
				else
					mm[r][c] = s + m[r][c];
	}
	///@}

	///@name Sized operator members
	///@{
	int num_rows() const
	{
		return m.num_rows();
	}
	int num_cols() const
	{
		return m.num_cols();
	}
	///@}
};

///@internal
///@brief Generic scalars object. Knows how to be added, knows how to deal with += and so on.
///See TooN::Ones
///@ingroup gInternal
template<class P> struct Operator<Internal::Scalars<P> >
{	
	///@name Scalable operator members
	///@{
	typedef P Precision;
	///@}

	const Precision s; ///<Value of the scalar being represented.  
	
	///@name Constructors
	///@{
	Operator(Precision s_)
		:s(s_){}

	Operator()
	{}
	///@}

	////////////////////////////////////////
	//
	// All applications for vector
	//
	///@name Operator members
	///@{

	template <int Size, typename P1, typename B1> 
	void eval(Vector<Size, P1, B1>& v) const
	{
		for(int i=0; i < v.size(); i++)
			v[i] = (P1)s;
	}

	template <int Size, typename P1, typename B1> 
	void plusequals(Vector<Size, P1, B1>& v) const
	{
		for(int i=0; i < v.size(); i++)
			v[i] += (P1)s;
	}

	template <int Size, typename P1, typename B1>
	void minusequals(Vector<Size, P1, B1>& v) const
	{
		for(int i=0; i < v.size(); ++i)
			v[i] -= (P1)s;
	}

	template <int Size, typename P1, typename B1> 
	Operator<Internal::ScalarsVector<Size,P1,B1,Precision> > add(const Vector<Size, P1, B1>& v) const
	{
		return Operator<Internal::ScalarsVector<Size,P1,B1,Precision> >(s, v, 0);
	}

	template <int Size, typename P1, typename B1> 
	Operator<Internal::ScalarsVector<Size,P1,B1,Precision> > rsubtract(const Vector<Size, P1, B1>& v) const
	{
		return Operator<Internal::ScalarsVector<Size,P1,B1,Precision> >(-s, v, 0);
	}

	template <int Size, typename P1, typename B1> 
	Operator<Internal::ScalarsVector<Size,P1,B1,Precision> > lsubtract(const Vector<Size, P1, B1>& v) const
	{
		return Operator<Internal::ScalarsVector<Size,P1,B1,Precision> >(s, v, 1);
	}

	////////////////////////////////////////
	//
	// All applications for matrix
	//

	template <int Rows, int Cols, typename P1, typename B1> 
	void eval(Matrix<Rows,Cols, P1, B1>& m) const
	{
		for(int r=0; r < m.num_rows(); r++)
			for(int c=0; c < m.num_cols(); c++)
				m[r][c] = s;
	}

	template <int Rows, int Cols, typename P1, typename B1> 
	void plusequals(Matrix<Rows,Cols, P1, B1>& m) const
	{
		for(int r=0; r < m.num_rows(); r++)
			for(int c=0; c < m.num_cols(); c++)
				m[r][c] += (P1)s;
	}

	template <int Rows, int Cols, typename P1, typename B1> 
	void minusequals(Matrix<Rows,Cols, P1, B1>& m) const
	{
		for(int r=0; r < m.num_rows(); r++)
			for(int c=0; c < m.num_cols(); c++)
				m[r][c] -= (P1)s;
	}

	template <int Rows, int Cols, typename P1, typename B1> 
	Operator<Internal::ScalarsMatrix<Rows,Cols,P1,B1,Precision> > add(const Matrix<Rows,Cols, P1, B1>& v) const
	{
		return Operator<Internal::ScalarsMatrix<Rows,Cols,P1,B1,Precision> >(s, v, 0);
	}


	template <int Rows, int Cols, typename P1, typename B1> 
	Operator<Internal::ScalarsMatrix<Rows,Cols,P1,B1,typename Internal::NegType<P>::Type> > rsubtract(const Matrix<Rows,Cols, P1, B1>& v) const
	{
		return Operator<Internal::ScalarsMatrix<Rows,Cols,P1,B1,typename Internal::NegType<P>::Type > >(-s, v, 0);
	}

	template <int Rows, int Cols, typename P1, typename B1> 
	Operator<Internal::ScalarsMatrix<Rows,Cols,P1,B1,Precision> > lsubtract(const Matrix<Rows,Cols, P1, B1>& v) const
	{
		return Operator<Internal::ScalarsMatrix<Rows,Cols,P1,B1,Precision> >(s, v, 1);
	}
	///@}
	////////////////////////////////////////
	//
	// Create sized versions for initialization
	//

	///@name Sizeable operators members
	///@{

	Operator<Internal::SizedScalars<Precision> > operator()(int size) const
	{
		return Operator<Internal::SizedScalars<Precision> > (s,size);
	}

	Operator<Internal::RCScalars<Precision> > operator()(int r, int c) const
	{
		return Operator<Internal::RCScalars<Precision> > (s,r,c);
	}
	///@}

	///@name Scalable operator members
	///@{
	template<class Pout, class Pmult> Operator<Internal::Scalars<Pout> > scale_me(const Pmult& m) const
	{
		return Operator<Internal::Scalars<Pout> >(s*m);
	}
	///@}
};

///@internal
///@brief Variant of the Operator<Internal::Scalars> object which holds a size to construct dynamic vectors or square matrices.
///@ingroup gInternal 
template<class P> struct Operator<Internal::SizedScalars<P> >: public Operator<Internal::Scalars<P> >
{
	using Operator<Internal::Scalars<P> >::s;
	///@name Sized operator members
	///@{
	const int my_size;
	int size() const {
		return my_size;
	}
	int num_rows() const {
		return my_size;
	}
	int num_cols() const {
		return my_size;
	}
	///@}

	///@name Constructors
	///@{
	Operator(P s, int sz)
		:Operator<Internal::Scalars<P> >(s),my_size(sz){}
	///@}
		
	///@name Scalable operator members
	///@{
	template<class Pout, class Pmult> Operator<Internal::SizedScalars<Pout> > scale_me(const Pmult& m) const
	{
		return Operator<Internal::SizedScalars<Pout> >(s*m, my_size);
	}
	///@}

private:
	void operator()(int);
	void operator()(int,int);
};

		
///@internal
///@brief Variant of Scalars (see TooN::Ones) which holds two sizes to construct dynamic matrices.
///@ingroup gInternal
template<class P> struct Operator<Internal::RCScalars<P> >: public Operator<Internal::Scalars<P> >
{
	using Operator<Internal::Scalars<P> >::s;

	///@name Operator members
	///@{
	const int my_rows, my_cols;
	int num_rows() const {
		return my_rows;
	}
	int num_cols() const {
		return my_cols;
	}
		
	Operator(P s, int r, int c)
		:Operator<Internal::Scalars<P> >(s),my_rows(r),my_cols(c)
	{}
		
	template<class Pout, class Pmult> Operator<Internal::RCScalars<Pout> > scale_me(const Pmult& m) const
	{
		return Operator<Internal::RCScalars<Pout> >(s*m, my_rows, my_cols);
	}

	///@}
private:
	void operator()(int);
	void operator()(int,int);
};


////////////////////////////////////////////////////////////////////////////////
//
// How to scale scalable operators
//
	
template<template<class> class Op, class Pl, class Pr> 
Operator<Op<typename Internal::MultiplyType<Pl, Pr>::type > >
operator*(const Pl& l, const Operator<Op<Pr> >& r)
{
	return r.template scale_me<typename Internal::MultiplyType<Pl, Pr>::type, Pl>(l); 
}

template<template<class> class Op, class Pl, class Pr> 
Operator<Op<typename Internal::MultiplyType<Pl, Pr>::type > >
operator*(const Operator<Op<Pl> >& l, const Pr&  r)
{
	return l.template scale_me<typename Internal::MultiplyType<Pl, Pr>::type>(r); 
}

template<template<class> class Op, class Pl, class Pr> 
Operator<Op<typename Internal::DivideType<Pl, Pr>::type > >
operator/(const Operator<Op<Pl> >& l, const Pr&  r)
{
	return l.template scale_me<typename Internal::MultiplyType<Pl, Pr>::type, Pl>(static_cast<typename Internal::DivideType<Pl,Pr>::type>(1)/r); 
}


template<class Op>
Operator<Op> operator-(const Operator<Op>& o)
{
	return o.template scale_me<typename Operator<Op>::Precision>(-1);
}

//Special case for negating One
template<template<class>class Op>
Operator<Op<DefaultPrecision> > operator-(const Operator<Op<Internal::One> >& o)
{
	return o.template scale_me<DefaultPrecision>(-1);
}

/**This function is used to add a scalar to every element of a vector or
   matrix. For example:
   @code
   Vector<3> v;
   ...
   ...
   v += Ones * 3; //Add 3 to every element of v;
   @endcode
   Both + and += are supported on vectors,matrices and slices.

   For construction of dynamic vectors and matrices, a size needs to be given:
   @code
       Vector<3> v_static = Ones;    
	   Vector<>  v_dynamic = Ones(3); //Construct a 3x1 vector full one 1s
       Matrix<3> m_static = Ones;     
	   Matrix<>  m_dynamic = Ones(3,4); //Construct a 3x4 matrix
   @endcode
   @ingroup gLinAlg
*/
static const Operator<Internal::Scalars<Internal::One> > Ones;


/**This function is used to initialize vectors and matrices to zero.
   For construction of dynamic vectors and matrices, a size needs to be given.
   For example:
   @code
       Vector<3> v_static = Zeros;
	   Vector<>  v_dynamic = Zeros(3); //Construct a 3x1 vector
       Matrix<3> m_static = Zeros;
	   Matrix<>  m_dynamic = Zeros(3,4); //Construct a 3x4 matrix
   @endcode
   @ingroup gLinAlg
*/
static Operator<Internal::Zero> Zeros;

/**This function is used to add a scalar to the diagonal of a matrix,
   or to construct matrices.
   For example:
   @code
   Matrix<3> v;
   ...
   ...
   Matrix<3> u = v  + Identity * 4;
   @endcode
   Both + and += are supported. For assignment, if the matrix is non-square,
   then all elements off the leading diagonal are set to zero.
   For construction of dynamic matrices, a size needs to be given:
   @code
       Matrix<3> m_static = Identity;     
	   Matrix<>  m_dynamic = Identity(3); //Construct a 3x3 matrix
   @endcode
   @ingroup gLinAlg
*/

static Operator<Internal::Identity<Internal::One> > Identity;

}