Logo Search packages:      
Sourcecode: libneuralnet version File versions  Download package

Matrix.hxx

//
// Matrix.hxx for  in
//
// Made by Samuel Pasquier
// Login   <pasqui_s@epita.fr>
//
// Started on  Sun May 19 16:45:14 2002 Samuel Pasquier
// Last update Tue May 28 23:41:22 2002 Samuel Pasquier
//

#ifndef _MATRIX_HXX_
# define _MATRIX_HXX_

# include <assert.h>
# include <math.h>
# include <exception/exception.hh>

//////////////////////////////////////////////////////////////////////
//
//                Constructors - Destructors
//
//////////////////////////////////////////////////////////////////////
template <class TypeName>
Matrix<TypeName>::Matrix(int lines, int cols)
{
  alloc_data(lines, cols);
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      _data[i][j] = 0;
}

template <class TypeName>
Matrix<TypeName>::Matrix(const Matrix<TypeName>& ref)
{
  alloc_data(ref.get_nb_lines(), ref.get_nb_cols());
  this->clone_same_size(ref);
}

template <class TypeName>
Matrix<TypeName>::Matrix(const std::vector<std::vector<TypeName> >& ref)
{
  assert(ref.size() > 0);
  alloc_data(ref.size(), ref[0].size());

  for (unsigned int i = 0; i < ref.size(); i++)
    for (unsigned int j = 0; j < ref[0].size(); j++)
      _data[i][j] = ref[i][j];
}

template <class TypeName>
Matrix<TypeName>::Matrix(const std::vector<TypeName>& ref)
{
  assert(ref.size() > 0);
  alloc_data(ref.size(), 1);

  for (unsigned int i = 0; i < ref.size(); i++)
    _data[i][0] = ref[i];
}

template <class TypeName>
Matrix<TypeName>::~Matrix()
{
  release_data();
}


template <class TypeName>
void Matrix<TypeName>::clone_same_size(const Matrix<TypeName>& ref)
{
  if ((_lines != ref.get_nb_lines()) || (_cols != ref.get_nb_cols()))
    throw Exception::Exception("Two Matrix have different size", HERE);
  if (!_data)
    throw Exception::Exception("Data are not malloc", HERE);
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      _data[i][j] = ref(i, j);
}

//////////////////////////////////////////////////////////////////////
//
//                Operator Overloading
//
//////////////////////////////////////////////////////////////////////
template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::operator -(const Matrix& rhs) const
{
  if ((_lines != rhs.get_nb_lines()) || (_cols != rhs.get_nb_cols()))
    throw Exception::Exception("Two Matrix have different size", HERE);
  if (!_data)
    throw Exception::Exception("Data are not malloc", HERE);
  Matrix<TypeName> tmp(_lines, _cols);
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      tmp.set_data(i, j, _data[i][j] - rhs(i, j));
  return tmp;
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::operator +(const Matrix& rhs) const
{
  if ((_lines != rhs.get_nb_lines()) || (_cols != rhs.get_nb_cols()))
    throw Exception::Exception("Two Matrix have different size", HERE);
  if (!_data)
    throw Exception::Exception("Data are not malloc", HERE);
  Matrix<TypeName> tmp(_lines, _cols);
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      tmp.set_data(i, j, _data[i][j] + rhs(i, j));
  return (tmp);
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::operator *(const Matrix& rhs) const
{
  if (_cols != rhs.get_nb_lines())
    throw Exception::Exception("Two Matrix have different size", HERE);
  if (!_data)
    throw Exception::Exception("Data are not malloc", HERE);
  Matrix<TypeName> tmp(_lines, rhs.get_nb_cols());
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < rhs.get_nb_cols(); j++)
      for (int k = 0; k < _cols; k++)
      tmp.set_data(i, j, tmp(i,j) + (_data[i][k] * rhs(k, j)));
  return (tmp);
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::operator *(TypeName t) const
{
  if (!_data)
    throw Exception::Exception("Data are not malloc", HERE);
  Matrix<TypeName> tmp(_lines, _cols);
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      tmp.set_data(i, j, _data[i][j] * t);
  return (tmp);
}

template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::operator +=(const Matrix& rhs)
{
  clone_same_size(*this + rhs);
  return *this;
}

template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::operator -=(const Matrix& rhs)
{
  clone_same_size(*this - rhs);
  return *this;
}

template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::operator *=(const Matrix& rhs)
{
  clone_same_size(*this * rhs);
  return *this;
}

template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::operator *=(TypeName t)
{
  clone_same_size(*this * t);
  return *this;
}

template <class TypeName>
Matrix<TypeName>&  Matrix<TypeName>::operator=(const Matrix<TypeName>& rhs)
{
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      _data[i][j] = rhs(i, j);

  return (*this);
}


//////////////////////////////////////////////////////////////////////
//
//                Data Accessors
//
//////////////////////////////////////////////////////////////////////
template <class TypeName>
const TypeName Matrix<TypeName>::operator ()(int x, int y) const
{
  if (!(((x >= 0) and (x < _lines)) and ((y >= 0) && (y <= _cols))))
    throw Exception::Exception("Matrix index out of bounds",
                         HERE);
  return (this->_data[x][y]);
}

template <class TypeName>
void Matrix<TypeName>::set_data(int x, int y, const TypeName& value)
{
  if (!(((x >= 0) and (x < _lines)) and ((y >= 0) && (y <= _cols))))
    throw Exception::Exception("Matrix index out of bounds",
                         HERE);
  this->_data[x][y] = value;
}

template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::insert(int x, int y, const Matrix<TypeName>& m)
{
  assert (x >= 0 and y >= 0 and
        x + m.get_nb_lines() <= _lines and
        y + m.get_nb_cols() <= _cols);
  for (int i = x; i < x + m.get_nb_lines(); i++)
    for (int j = y; j < y + m.get_nb_cols(); j++)
      _data[i][j] = m(i - x, j - y);
  return *this;
}

template <class TypeName>
int Matrix<TypeName>::get_nb_lines() const
{
  return _lines;
}

template <class TypeName>
int Matrix<TypeName>::get_nb_cols() const
{
  return _cols;
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::get_line(int line) const
{
  if (line > _lines)
    throw Exception::Exception("Matrix index out of bounds", HERE);
  Matrix<TypeName> c(1, _cols);
  for (int i = 0; i < _cols; i++)
    c.set_data(0, i, _data[line][i]);
  return c;
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::get_col(int col) const
{
  if (col > _cols)
    throw Exception::Exception("Matrix index out of bounds", HERE);
  Matrix<TypeName> c(_lines, 1);
  for (int i = 0; i < _lines; i++)
    c.set_data(i, 0, _data[i][col]);
  return c;
}

template <class TypeName>
Matrix<TypeName>
Matrix<TypeName>::get_cols(int nb_cols) const
{
  Matrix<TypeName> result(_lines, nb_cols);

  for (int j = 0; j < nb_cols; j++)
    for (int i = 0; i < _lines; i++)
      result.set_data(i, j, _data[i][j]);
  return result;
}

template <class TypeName>
Matrix<TypeName>
Matrix<TypeName>::get_lines(int nb_lines) const
{
  Matrix<TypeName> result(nb_lines, _cols);

  for (int j = 0; j < _cols; j++)
    for (int i = 0; i < nb_lines; i++)
      result.set_data(i, j, _data[i][j]);
  return result;
}


//////////////////////////////////////////////////////////////////////
//
//                Matrix operations
//
//////////////////////////////////////////////////////////////////////


template <class TypeName>
TypeName Matrix<TypeName>::trace() const
{
  TypeName res = 0;

  assert(_lines == _cols);
  for (int i = 0; i < _lines; i++)
    res += _data[i][i];
  return (res);
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::abs() const
{
  Matrix<TypeName> res(_lines, _cols);

  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      if (_data[i][j] < 0)
      res.set_data(i, j, -_data[i][j]);
      else
      res.set_data(i, j, _data[i][j]);
  return res;
}


template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::transpose() const
{
  int t_lines = _cols;
  int t_cols = _lines;
  Matrix<TypeName> t(t_lines, t_cols);

  for (int i = 0; i < t_lines; i++)
    for (int j = 0; j < t_cols; j++)
      t.set_data(i, j, _data[j][i]);
  return t;
}

template <class TypeName>
Matrix<double> Matrix<TypeName>::inverse() const
{
  Matrix<TypeName> ak = get_col(0);
  Matrix<double> Aprime(_cols, _lines);
  if (ak.norme() > 1e-5)
    Aprime.insert(0, 0, ak.transpose() * (1 / (double)ak.norme()));

  for (int k = 1; k < _cols; k++)
    {
      ak = get_col(k);
      Matrix<double> d = Aprime.get_lines(k) * ak;
      Matrix<double> c = ak - get_cols(k) * d;
      Matrix<double> b(1, _lines);
      if (c.norme() < 1e-5)
      b.insert(0, 0, ((Aprime.get_lines(k)).transpose() * d * (1 / (1 + d.norme()))).transpose());
      else
      b.insert(0, 0, c.transpose() * (1 / (double)c.norme()));
      Aprime.insert(0, 0, (Aprime.get_lines(k) - d * b));
      Aprime.insert(k, 0, b);
    }
  return Aprime;
}

template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::identity()
{
  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      if (i == j)
      _data[i][j] = 1;
      else
      _data[i][j] = 0;
  return *this;
}

/*
** renvoi un vecteur des moyennes de chaque colonne (util pour obtenir
** la matrice reduite)
*/
template <class TypeName>
Matrix<float> Matrix<TypeName>::center()
{
  Matrix<float> moy(_cols, 1);

  for (int j = 0; j < _cols; j++)
    {
      TypeName tmp = 0;
      for (int i = 0; i < _lines; i++)
      tmp += _data[i][j];
      float cur_moy = tmp / (float)_lines;
      moy.set_data(j, 0, cur_moy);
      for (int i = 0; i < _lines; i++)
      _data[i][j] -= cur_moy;
    }
  return moy;
}

/*
** renvoi un vecteur des ecart-type de chaque colonne
*/
template <class TypeName>
Matrix<float> Matrix<TypeName>::reduce(Matrix<float>& moy)
{
  Matrix<float> et(_cols, 1);

  for (int j = 0; j < _cols; j++)
    {
      TypeName tmp = 0;
      for (int i = 0; i < _lines; i++)
      tmp += (_data[i][j] - moy(j, 0)) * (_data[i][j] - moy(j, 0));
      tmp /= _lines;
      et.set_data(j, 0, sqrt(tmp));
      for (int i = 0; i < _lines; i++)
      if (et(j, 0))
        _data[i][j] /= et(j, 0);
    }
  return et;
}

template <class TypeName>
Matrix<float> Matrix<TypeName>::get_varcov() const
{
  return (this->transpose() * (*this) * (1 / (float)_lines));
}

template <class TypeName>
float Matrix<TypeName>::trace_of_varcov() const
{
  float result = 0;

  for (int i = 0; i < _lines; i++)
    for (int j = 0; j < _cols; j++)
      result += _data[i][j] * _data[i][j];
  return result / (float)_lines;
}

//This function assumes the _data attribute are already
//deleted (with release_data)
template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::alloc_data(int lines, int cols)
{
  _data = new TypeName * [lines];
  for (int i = 0; i < lines; i++)
    _data[i] = new TypeName[cols];
  _lines = lines;
  _cols = cols;
  return *this;
}

template <class TypeName>
Matrix<TypeName>& Matrix<TypeName>::release_data()
{
  assert(_data);
  for (int i = 0; i < _lines; i++)
    delete [] _data[i];
  delete [] _data;
  _data = 0;
  _lines = _cols = 0;
  return *this;
}


//////////////////////////////////////////////////////////////////////
//
//                Streamable capability
//
//////////////////////////////////////////////////////////////////////
template <class TypeName>
std::ostream& Matrix<TypeName>::print(std::ostream& os) const
{
  os << std::endl;
  for (int i = 0; i < _lines; i++)
    {
      for (int j = 0; j < _cols; j++)
      os << "[ " << _data[i][j] << " ]";
      os << std::endl;
    }
  return (os);
}

template <class TypeName>
std::ostream& operator <<(std::ostream& os, const Matrix<TypeName> & m)
{
  return (m.print(os));
}

template <class TypeName>
TypeName Matrix<TypeName>::norme() const
{
  if (_cols != 1)
    throw Exception::Exception("Two Matrix have different size", HERE);
  if (!_data)
    throw Exception::Exception("Data are not malloc", HERE);
  TypeName res = 0;
  for (int i = 0; i < _lines; i++)
    res += _data[i][0] * _data[i][0];
  return (res);
}

template <class TypeName>
TypeName Matrix<TypeName>::norme(const Matrix& rhs) const
{
  if ((_lines != rhs.get_nb_lines()) || (_cols != 1))
    throw Exception::Exception("Two Matrix have different size", HERE);
  if (!_data)
    throw Exception::Exception("Data are not malloc", HERE);
  TypeName res = 0;
  for (int i = 0; i < _lines; i++)
    res += _data[i][0] * rhs(i, 0);
  return (res);
}

template <class TypeName>
void Matrix<TypeName>::prepare_next_eigen_vector(Matrix<TypeName>& current_q,
                                     double current_lambda)
{
  (*this) -= current_q * current_q.transpose() * current_lambda;
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::get_max_eigen_vector() const
{
  Matrix<TypeName> result(_lines, 1);
  result.set_data(1, 0, 0);
  for (int i = 1; i < result.get_nb_lines(); i++)
    result.set_data(i, 0, 1);
  TypeName l , l1;
  l = 0.0;
  Matrix<TypeName> y(_lines, 1);
  do
    {
      l1 = l;
      y.clone_same_size((*this) * result);
      l = sqrt(y.norme());
      for (int i = 0; i < _lines; i++)
      result.set_data(i, 0, y(i, 0) / l);
    }
  while (fabs(double(l - l1)) > 1e-10);
  return result;
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::get_eigen_value(Matrix<TypeName>& eigenvector) const
{
  return eigenvector.transpose() * (*this) * eigenvector;
}

template <class TypeName>
Matrix<TypeName> Matrix<TypeName>::get_max_eigen_vector_deflation() const
{
  Matrix<TypeName> v(_lines, 1);
  Matrix<TypeName> w(_lines, 1);

  for (int i = 0; i < v.get_nb_lines(); i++)
    {
      v.set_data(i, 0, random() % 3 + 1);
      w.set_data(i, 0, 1);
    }
  TypeName l , l1, n;
  n = v.norme();
  for (int i = 0; i < v.get_nb_lines(); i++)
    v.set_data(i, 0, v(i, 0) / n);
  l = 0.0;
  Matrix<TypeName> y(_lines, 1);
  std::cout << "v =" << v << std::endl;
  do
    {
      l1 = l;
      y.clone_same_size((*this) * w);
      //      std::cout << "y = " << y << std::endl;
      l = y.norme(v);
      std::cout << "l = " << l << std::endl;
      for (int i = 0; i < _lines; i++)
            w.set_data(i, 0, y(i, 0) / l);
      std::cout << "w = " << w << std::endl
            << "l - l1 = " << l - l1 << std::endl;

    }
  while (fabs(double(l - l1)) > 1e-6);
  return w;
}


#endif /* !_MATRIX_HXX_ */



Generated by  Doxygen 1.6.0   Back to index