/* -*- mona-c++ -*-
* Copyright (c) Leipzig, Madrid 2004 - 2008
* Max-Planck-Institute for Human Cognitive and Brain Science
* Max-Planck-Institute for Evolutionary Anthropology
* BIT, ETSI Telecomunicacion, UPM
*
* 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 2 of the License, or
* (at your option) any later version.
*
* This program 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 program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
// $Id: Vector.hh 740 2005-11-22 12:57:21Z write1 $
/*! \brief Implements n-dimensional vectors, nxm - matrices, and some special squared matrices
This provides a n-Dim-vector class, which uses reference-counting
Datafield to avoid unneccesary copying of the data-fields
\todo make the storage a policy
\file Vector.hh
\author Gert Wollny <wollny@cbs.mpg.de>
*/
#ifndef __MONA_VECTOR_HH
#define __MONA_VECTOR_HH 1
#if HAVE_CONFIG_H
#include <config.h>
#endif
#include <cstring>
#include <stdexcept>
#include <cassert>
#include <cmath>
#include <cstdio>
#include <numeric>
// MONA specific
#include <libmona/monaTypes.hh>
#include <libmona/RefData.hh>
#include <libmona/3DVector.hh>
#include <libmona/3DAlgorithm.hh>
#ifndef TRACE
#define TRACE(a,b)
#endif
namespace mona {
template <class T> class TVector {
typedef TRefData<T> RefData;
protected:
TVector(RefData *data);
void release_data();
public:
typedef T Element;
TVector();
/** constructor to create a vector of dize _dim and iniztialize its elemtns to zero.
\param _dim size of vector
*/
TVector(int _dim);
/** copy constructor, takes onyl a reference to the data (may change to policy)
\param org the vector to be copied
*/
TVector(const TVector& org);
/** virtual destructor
\remark why is this virtual?
*/
virtual ~TVector();
/** Assignment operator takes only the reference, may chaneg with policy
\param org TVector to assign from
*/
TVector& operator = (const TVector& org);
/** Euclidian norm of vector */
double norm()const;
/** Square of Euclidian norm of vector */
double norm2()const;
/** Set all elements to T() (may change with policy */
virtual void clear();
/** read only index operator */
const Element& operator [] (unsigned int i) const;
/** read/write index operator */
Element& operator [] (unsigned int i);
/** in-place addition
\param v the vector do be added
\todo change return value
*/
void operator += (const TVector<T>& v);
/** in-place subtraction
\param v the vector do be subtracted
\todo change return value
*/
void operator -= (const TVector<T>& v);
/** in-place multiplication with scalar
\param d the scalar to multiply by
\todo change return value
*/
void operator *= (float d);
/** in-place division by scalar
\param d the scalar to divide by
\todo change return value
*/
void operator /= (float d);
TVector<T> operator + (const TVector<T>& v)const;
TVector<T> operator - (const TVector<T>& v)const;
double operator * (const TVector<T>& v)const;
TVector<T> operator * (float f)const;
TVector<T> operator / (float f)const;
/// Ensures, that \a this uses a single referenced data block
virtual void make_single_ref(){check_ref_count();};
bool operator == (const TVector& v)const;
bool operator != (const TVector& v)const;
/// returns dimension of vector
unsigned int size() const;
private:
void check_ref_count(); // if other objects share this data, make a copy before changing
const static Element Zero;
RefData *_M_data;
};
/** A Matrix definrd as vector of vectors of type T.
\todo decide if it is better not to derive from \a TVector but store data in an TRefData policy
*/
template <class V> class TMatrix:public TVector< TVector <V> >{
public:
typedef TVector<V> Vector;
typedef V Element;
/** Constructor
\param iRows_ number of rows in mtrix
\param iColums_ number of columns in matzrix
*/
TMatrix(int iRows_, int iColums_);
/// returns number of rows
unsigned int get_rows() const;
/// returns number of columns
unsigned int get_cols() const;
/// Ensures, that \a this uses a single referenced data block
virtual void make_single_ref();
/// makes the matrix a Zero - Matrix
virtual void clear();
/** get a column vector
\param i requested column
\returns TVector of column i
*/
Vector get_column(unsigned int i) const;
/** read-write element access operator
\param i requested row
\param j requested column
\returns reference to matrix element at (i,j)
*/
Element& operator()(unsigned int j,unsigned int i);
/** read-only element access operator
\param i requested row
\param j requested column
\returns const reference to matrix element at (i,j)
*/
const Element& operator()(unsigned int j,unsigned int i)const;
/** returns the transposed of this matrix
\todo defined this outside of the class
*/
virtual TMatrix<V> transpose() const;
/** set a column of the matrix
\param i number of column to be set
\param c column vector to set to
*/
void set_column(unsigned int i, const Vector& c);
/** Matrix-multiplication operator
\todo move this outside of class
*/
TMatrix<V> operator * (const TMatrix<V>& m);
/** multiplication of matrix with vector
\todo move this outside of class
*/
Vector operator * (const Vector& a)const;
};
template <class T> class TQuadMatrix: public TMatrix<T> {
public:
typedef TVector<T> Vector;
typedef T Element;
/// a square matrix
TQuadMatrix(int Dim);
/// returns always \a false
virtual bool is_decomp()const;
/// not implemented, a \a TQuadMatrix can not be set to decomposed per se
virtual void set_decomp();
};
template <class T> class TLRMatrix: public TQuadMatrix<T> {
TVector<int> Perm;
public:
typedef TVector<T> Vector;
typedef T Element;
/**
Create a LR-decomposed Matrix from an \a TQuadMatrix
*/
TLRMatrix(const TQuadMatrix<T>& org);
int get_perm(int k)const;
virtual bool is_decomp()const;
Vector operator *(const Vector& v) const;
};
template <class T>
class TSymmMatrix: public TQuadMatrix<T> {
bool bIsDecomp;
public:
typedef TVector<T> Vector;
typedef T Element;
TSymmMatrix(int Dim);
// virtual bool invert();
Element& operator()(unsigned int j,unsigned int i);
const Element& operator()(unsigned int j,unsigned int i)const;
virtual TMatrix<T> Transpose();
virtual bool is_decomp()const;
virtual void set_decomp();
Vector operator *(const Vector& v) const;
};
typedef TMatrix<float> CFMatrix;
typedef TVector<float> CFVector;
typedef TVector<C3DFVector> CSolverVector;
typedef TQuadMatrix<float> CFQuadMatrix;
typedef TLRMatrix<float> CFLRMatrix;
typedef TSymmMatrix<float> CFSymmMatrix;
//typedef TRefData<float> CFData;
//typedef TRefData<C3DFVector> CSolverData;
template <class V>
unsigned int TMatrix<V>::get_rows() const
{
return this->size();
}
template <class V>
unsigned int TMatrix<V>::get_cols() const
{
return (*this)[0].size();
}
template <class V>
void TMatrix<V>::make_single_ref()
{
for (size_t i = 0; i < this->size(); i++) {
(*this)[i].make_single_ref();
}
TVector< TVector <V> >::make_single_ref();
}
template <class V>
void TMatrix<V>::clear()
{
for (size_t i = 0; i < this->size(); i++) {
(*this)[i].clear();
}
}
template <class V>
typename TMatrix<V>::Element&
TMatrix<V>::operator()(unsigned int j,unsigned int i)
{
return (*this)[j][i];
}
template <class V>
const typename TMatrix<V>::Element&
TMatrix<V>::operator()(unsigned int j,unsigned int i)const
{
return (*this)[j][i];
}
template <class T>
bool TQuadMatrix<T>::is_decomp()const
{
return false;
}
template <class T>
bool TLRMatrix<T>::is_decomp()const
{
return true;
}
template <class T>
int TLRMatrix<T>::get_perm(int k)const
{
return Perm[k];
}
template <class T>
typename TSymmMatrix<T>::Element&
TSymmMatrix<T>::operator()(unsigned int j,unsigned int i)
{
if (i<=j) {
return (*this)[i][j];
}else{
return (*this)[j][i];
}
}
template <class T>
const typename TSymmMatrix<T>::Element&
TSymmMatrix<T>::operator()(unsigned int j,unsigned int i)const
{
if (i<=j) {
return (*this)[i][j];
}else{
return (*this)[j][i];
}
}
template <class T>
TMatrix<T> TSymmMatrix<T>::Transpose()
{
fprintf(stderr,"WARNING:TSymmMatrix<T>::Transpose() returnes *this, but only a TMatrix\n");
return *this;
}
template <class T>
bool TSymmMatrix<T>::is_decomp()const
{
return bIsDecomp;
}
template <class T>
void TSymmMatrix<T>::set_decomp()
{
bIsDecomp = true;
}
template <class T> void TVector<T>::check_ref_count()
{ // if other objects share this data, make a copy before changing
if (_M_data->get_refcount() > 1) {
_M_data->release();
_M_data = new RefData(*_M_data);
_M_data->add_ref();
}
}
template <typename T>
TVector<T>::TVector(RefData *data):_M_data(data)
{
}
template <class T> TVector<T>::TVector():_M_data(NULL)
{
}
template <class T> TVector<T>::TVector(int _dim)
{
_M_data = new RefData(_dim);
_M_data->add_ref();
}
template <class T> TVector<T>::TVector(const TVector& org)
{
_M_data = org._M_data;
if (_M_data) {
_M_data->add_ref();
}
}
template <class T>
void TVector<T>::clear()
{
check_ref_count();
if (_M_data) {
_M_data->clear();
}
}
template <typename T>
void TVector<T>::release_data()
{
_M_data = NULL;
}
template <class T>
TVector<T>::~TVector()
{
if (_M_data)
_M_data->release();
}
template <class T>
TVector<T>& TVector<T>::operator = (const TVector& org)
{
if (&org == this)
return *this;
if (_M_data)
_M_data->release();
_M_data = org._M_data;
if (_M_data) {
_M_data->add_ref();
}
return *this;
}
template <class T>
const typename TVector<T>::Element&
TVector<T>::operator[] (unsigned int i) const
{
return _M_data->at(i);
}
template <class T>
typename TVector<T>::Element&
TVector<T>::operator [] (unsigned int i) {
check_ref_count();
assert(i < size());
return (*_M_data)[i];
}
template <class T> void TVector<T>::operator += (const TVector& v)
{
assert(v.size() == size());
check_ref_count();
typename TRefData<T>::iterator ime = _M_data->begin();
typename TRefData<T>::iterator eme = _M_data->end();
typename TRefData<T>::iterator vi = v._M_data->begin();
while (ime != eme) {
*ime += *vi;
++ime;
++vi;
}
}
template <class T> void TVector<T>::operator -= (const TVector& v)
{
assert(v.size() == size());
check_ref_count();
typename TRefData<T>::iterator ime = _M_data->begin();
typename TRefData<T>::iterator eme = _M_data->end();
typename TRefData<T>::iterator vi = v._M_data->begin();
while (ime != eme) {
*ime -= *vi;
++ime;
++vi;
}
}
template <class T> bool TVector<T>::operator == (const TVector& v)const
{
assert(v.size() == size());
return *_M_data == *v._M_data;
}
template <class T> bool TVector<T>::operator != (const TVector& v)const
{
return !(*this == v);
}
template <class T> void TVector<T>::operator *= (float d)
{
check_ref_count();
typename TRefData<T>::iterator ime = _M_data->begin();
typename TRefData<T>::iterator eme = _M_data->end();
while (ime != eme) {
*ime *= d;
++ime;
}
}
template <class T> void TVector<T>::operator /= (float d)
{
assert(d);
d = 1/d;
*this *= d;
}
template <class T> TVector<T> TVector<T>::operator + (const TVector<T>& v)const
{
assert(v.size() == size());
TVector <T> retVal(*this);
retVal += v;
return retVal;
}
template <class T> TVector<T> TVector<T>::operator - (const TVector<T>& v)const
{
assert(v.size() == size());
TVector <T> retVal(*this);
retVal -= v;
return retVal;
}
template <class T> double TVector<T>::operator * (const TVector<T>& v)const
{
double result = 0.0;
assert(v.size() == size());
typename TRefData<T>::iterator ime = _M_data->begin();
typename TRefData<T>::iterator eme = _M_data->end();
typename TRefData<T>::iterator vi = v._M_data->begin();
while (ime != eme) {
result += *ime * *vi;
++ime;
++vi;
}
return result;
}
template <class T> TVector<T> TVector<T>::operator * (float e) const
{
TVector <T> retVal(size());
size_t i = 0;
while (i < _M_data->size()){
retVal[i] = _M_data->at(i) * e;
i++;
}
return retVal;
}
template <class T> TVector<T> TVector<T>::operator / (float f)const
{
assert(f);
f = 1/f;
return *this * f;
}
template <class T> double TVector<T>::norm()const
{
return sqrt(*this * *this);
}
template <class T> unsigned int TVector<T>::size() const
{
return _M_data->size();
}
template <class T> double TVector<T>::norm2()const
{
return *this * *this;
}
template <class T>
const typename TVector<T>::Element
TVector<T>::Zero = TVector<T>::Element();
template <class V>
TMatrix<V>::TMatrix(int iRows, int iCols):TVector< TVector<V> >(iRows)
{
for (int i = 0; i < iRows; i++) {
(*this)[i] = Vector(iCols);
}
}
template <class V>
typename TMatrix<V>::Vector
TMatrix<V>::operator * (const Vector& v)const
{
assert(get_cols() == v.size());
Vector Result(this->size());
for (size_t i = 0; i < this->size(); i++) {
Result[i] = (*this)[i] * v;
}
return Result;
}
template <class V>
typename TMatrix<V>::Vector
TMatrix<V>::get_column(unsigned int i) const
{
assert(i < get_cols());
Vector Result(get_rows());
for (int j = 0; j < get_rows(); j++) {
Result[j] = (*this)[j][i];
}
return Result;
}
template <class V>
void TMatrix<V>::set_column(unsigned int i, const Vector& c)
{
assert(i < get_cols());
assert(c.size() == get_rows());
for (int j = 0; j < get_rows(); j++) {
(*this)[j][i] = c[j];
}
}
template <class V>
TMatrix<V> TMatrix<V>::transpose()const
{
TMatrix<V> Result(get_cols(),get_rows());
for (size_t i = 0; i< get_cols(); i++){
for (size_t j = 0;j< get_rows(); j++){
Result(i,j) = (*this)(j,i);
}
}
return Result;
}
template <class V>
TMatrix<V> TMatrix<V>::operator * (const TMatrix<V>& m)
{
assert(m.get_rows() == get_cols() );
TMatrix<V> result(get_rows(),m.get_cols());
for (int i = 0; i < get_cols(); i++) {
result.set_column(i,(*this) * m.get_column(i));
}
return result;
}
template <class T>
TQuadMatrix<T>::TQuadMatrix(int Dim):TMatrix<T>(Dim,Dim)
{
}
template <class T>
TLRMatrix<T>::TLRMatrix(const TQuadMatrix<T>& org):
TQuadMatrix<T>(org),Perm(org.size())
{
Vector d(this->size());
if ( org.is_decomp() ) {
throw std::range_error("already decomposited type");
}
for (int k = 0; k < this->size(); k ++) {
Perm[k] = k;
float zmax = 0.0;
for (int j = 0; j < this->size(); j++) {
zmax = Max(zmax, float(fabs((*this)(k,j))));
}
if (zmax == 0.0) {
throw std::range_error("Matrix singular");
}
d[k] = zmax;
}
// double signdet = 1.0;
for (int k = 0; k < this->size()-1; k++) {
double piv = fabs((*this)(k,k)/ d[k]);
int j0 = k;
for (int j = k+1; j < this->size(); j++) {
double tmp = fabs((*this)(j,k)/d[j]);
if (piv < tmp) {
piv = tmp; j0 = j;
}
}
if (piv < MACH_EPS) {
throw std::range_error("matrix numeric singular");
}
if (j0 != k) {
do_swap(&Perm[j0],&Perm[k]);
do_swap(&d[j0],&d[k]);
Vector hlp = (*this)[j0];
(*this)[j0] = (*this)[k] ;
(*this)[k] = hlp;
}
for (int j = k+1; j < this->size(); j++) {
if ((*this)(j,k) != 0.0){
(*this)(j,k) /= (*this)(k,k);
Element tmp = (*this)(j,k);
for (int m = k+1; m < this->size(); m++){
(*this)(j,m) -= tmp * (*this)(k,m);
}
}
}
}// end k;
if (fabs( (*this)(this->size()-1, this->size()-1)) < MACH_EPS) {
throw std::range_error("matrix numeric singular");
}
}
template <class Vector,class T>
Vector operator / (const Vector& b, const TLRMatrix<T>& lu)
{
// More exactly this is solving LR * x = b
Vector x(b.size());
for (int k = 0; k < lu.size(); k++) {
x[k] = b[lu.get_perm(k)];
for (int j = 0; j < k; j++) {
x[k] -= lu(k,j) * x[j];
}
}
for (int k = lu.size()-1; k >= 0 ; k--) {
typename TLRMatrix<T>::Element sum = 0;
for (int j = k+1; j < lu.size(); j++) {
sum += lu(k,j) * x[j];
}
x[k] = (x[k] - sum)/lu(k,k);
}
return x;
}
template <class T>
typename TLRMatrix<T>::Vector
TLRMatrix<T>::operator *(const Vector& v) const
{
// we are not interested in the multiplication of a LR
// decomposit with a vector
throw std::range_error("operator not implemented");
}
template <class T>
TSymmMatrix<T>::TSymmMatrix(int dim):
TQuadMatrix<T>(dim),bIsDecomp(false)
{
}
template <class T>
typename TSymmMatrix<T>::Vector
TSymmMatrix<T>::operator *(const Vector& v) const
{
Vector Result(this->size());
if (bIsDecomp) {
Vector Help(this->size());
//Help = Rt * v
for (int i = 0; i < this->size(); i++) {
typename Vector::Element help = Vector::Element();
for (int j = i; j < this->size(); j++){
help += (*this)(i,j) * v[j];
}
Help[i] = help;
}
// Result = R * Help
for (int i = 0; i < this->size(); i++) {
typename Vector::Element help = Vector::Element();
for (int j = 0; j <=i; j++){
help += (*this)(i,j) * Help[j];
}
Result[i] = help;
}
}else{
for (int i = 0; i < this->size(); i++) {
typename Vector::Element help = Vector::Element();
for (int j = 0; j < this->size(); j++){
help += (*this)(i,j) * v[j];
}
Result[i] = help;
}
}
return Result;
}
} // namespace mona
#endif
/* CVS LOG
$Log$
Revision 1.7 2005/08/18 12:24:04 tittge
bug in line 647 removed
Revision 1.6 2005/06/29 13:22:20 wollny
switch to version 0.7
Revision 1.2 2005/05/12 08:13:42 gerddie
clean up the configuration and remove FLOAT
Revision 1.1.1.1 2005/03/17 13:44:20 gerddie
initial import
Revision 1.5 2004/08/25 09:08:31 wollny
added an emacs style comment to all source files
Revision 1.4 2004/07/13 09:16:16 wollny
g++ 3.4 compile
Revision 1.3 2004/06/03 09:57:32 wollny
Changed (hopefully) all instancable class names to Cxxxxx
Revision 1.2 2004/02/12 14:00:25 tittge
Major addaptions from libmia0
Revision 1.1.1.1 2004/02/11 18:18:05 tittge
start project
*/