NICE_VisLearning/math/algebra/GMSparse2.cpp

#include "GMSparse2.h"
#include <assert.h>

using namespace OBJREC;
using namespace NICE;
using namespace std;

GMSparse2::GMSparse2 ( uint _rows, uint _cols ) : m_rows(_rows), m_cols(_cols) 
{
	resize(_rows, _cols);
}

GMSparse2::GMSparse2 ( const NICE::Matrix & iceA, double epsilon)
{
	m_rows = iceA.rows();
	m_cols = iceA.cols();
	newvectors = true;
	for(uint r = 0; r < m_rows; r++)
	{
		SparseVector *tmp = new SparseVector(m_cols);
		for(uint c = 0; c < m_cols; c++)
		{
			if(fabs(iceA(r ,c)) > epsilon)
				(*tmp)[c] = iceA(r, c);
		}
		A.push_back(tmp);
	}
}

GMSparse2::~GMSparse2()
{
	clear();
}

SparseVector & GMSparse2::operator[](int i) const
{
	return *A[i];
}

void GMSparse2::addRow(const NICE::Vector &x)
{
	SparseVector *v = new SparseVector(x);
	m_cols = x.size();
	v->setDim(m_cols);
	A.push_back(v);
	++m_rows;
}
		
void GMSparse2::addRow(SparseVector *x)
{
	newvectors = false;
	A.push_back(x);
	++m_rows;
	m_cols = x->getDim();
}
		
void GMSparse2::clear()
{
	if(newvectors)
	{
		for(uint i = 0; i < A.size(); i++)
		{
			delete A[i];
		}
	}
	A.clear();
}
		
void GMSparse2::resize( int _rows, int _cols)
{
	m_rows = _rows;
	m_cols = _cols;
	newvectors = true;
	for(uint i = 0; i < m_rows; i++)
	{
		SparseVector *tmp = new SparseVector(m_cols);
		A.push_back(tmp);
	}
}
		
void GMSparse2::mult2(GMSparse2 &y, GMSparse2 &out, bool transpx, bool transpy)
{
	int rowsx, colsx, rowsy, colsy;
	if(transpx)
	{
		colsx = m_rows;
		rowsx = m_cols;
	}
	else
	{
		rowsx = m_rows;
		colsx = m_cols;
	}
	
	if(transpy)
	{
		colsy = y.rows();
		rowsy = y.cols();
	}
	else
	{
		rowsy = y.rows();
		colsy = y.cols();
	}
	
	if(rowsy != colsx)
	{
		cout << "GMSparse::mult matrix sizes missmatched" << endl;
		exit(1);
	}
	
	out.resize(rowsx,colsy);

	double xval, yval, val, val2;
	
	for(int c = 0; c < colsy; c++)
	{
		for(int i = 0; i < rowsy; i++)
		{
			if(transpy)
				yval = y[c].get(i);
			else
				yval = y[i].get(c);
			
			if(yval != 0.0)
			{
				for(int r = 0; r < rowsx; r++)
				{
					if(transpx)
						xval = (*A[i]).get(r);
					else
						xval = (*A[r]).get(i);
					
					val = out[r].get(c);
					val2 = val+xval*yval;
					if(fabs(val2) >  10e-10)
						out[r][c] = val2;
					else if(val != 0.0)
						out[r].erase(c);
				}
			}
		}
	}
}
		
void GMSparse2::mult(GMSparse2 &y, GMSparse2 &out, bool transpx, bool transpy)
{
	int rowsx, colsx, rowsy, colsy;
	if(transpx)
	{
		colsx = m_rows;
		rowsx = m_cols;
	}
	else
	{
		rowsx = m_rows;
		colsx = m_cols;
	}
	
	if(transpy)
	{
		colsy = y.rows();
		rowsy = y.cols();
	}
	else
	{
		rowsy = y.rows();
		colsy = y.cols();
	}
	
	if(rowsy != colsx)
	{
		cout << "GMSparse::mult matrix sizes missmatched" << endl;
		exit(1);
	}
	
	out.resize(rowsx,colsy);

	double xval, yval, val;
	for(int r = 0; r < rowsx; r++)
	{
	
		for(int c = 0; c < colsy; c++)
		{
			val = 0.0;
			for(int i = 0; i < rowsy; i++)
			{
				if(transpx)
					xval = (*A[i]).get(r);
				else
					xval = (*A[r]).get(i);
				
				if(transpy)
					yval = y[c].get(i);
				else
					yval = y[i].get(c);
				
				val+=xval * yval;
			}

			if(fabs(val) >  10e-7)
				out[r][c] = val;
			
			if(!finite(val*val))
			{
				cout << "val " << val << endl;
				for(int i = 0; i < rowsy; i++)
				{
					if(transpx)
						xval = (*A[i]).get(r);
					else
						xval = (*A[r]).get(i);
				
					if(transpy)
						yval = y[c].get(i);
					else
						yval = y[i].get(c);
				
					val+=xval * yval;
					cout << " xval: " << xval << " yval " << yval << endl;
				}
				getchar();
			}
		}
	}
	
	
}

void GMSparse2::mult(SparseVector &y, GMSparse2 &out, bool transpx, bool transpy)
{
	int rowsx, colsx, rowsy, colsy;
	if(transpx)
	{
		colsx = m_rows;
		rowsx = m_cols;
	}
	else
	{
		rowsx = m_rows;
		colsx = m_cols;
	}
	
	if(transpy)
	{
		colsy = 1;
		rowsy = y.getDim();
	}
	else
	{
		rowsy = 1;
		colsy = y.getDim();
	}
	
	if(rowsy != colsx)
	{
		cout << "GMSparse::mult matrix sizes missmatched" << endl;
		exit(1);
	}
	
	out.resize(rowsx,colsy);

	double xval, yval, val;
	for(int r = 0; r < rowsx; r++)
	{
	
		for(int c = 0; c < colsy; c++)
		{
			val = 0.0;
			for(int i = 0; i < rowsy; i++)
			{
				if(transpx)
					xval = (*A[i]).get(r);
				else
					xval = (*A[r]).get(i);
				
				if(transpy)
					yval = y.get(i);
				else
					yval = y.get(c);
				
				val += xval * yval;
			}

			if(fabs(val) >  10e-7)
				out[r][c] = val;
		}
	}
}

void GMSparse2::multiply ( NICE::Vector & y, const NICE::Vector & x ) const
{
	//FIXME: einbauen
}

void GMSparse2::store (ostream & os, int format) const
{
	os << m_rows << " " << m_cols << " " << A.size() << endl;
	for(int i = 0; i < (int) A.size(); i++)
	{
		A[i]->store(os, format);
	}
}

void GMSparse2::restore (istream & is, int format)
{
	A.clear();
	is >> m_rows;
	is >> m_cols;
	int size;
	is >> size;
	for(int i = 0; i < size; i++)
	{
		SparseVector *s = new SparseVector(m_cols);
		s->restore(is, format);
		A.push_back(s);
	}
}

void GMSparse2::setDel(bool del)
{
	newvectors = del;
}