// BigMatrix.h: interface for the CBigMatrix class.
//
// The big matrix is a variable-size matrix.
// THe main purpose of this class is facilitation of the problem of
// solving a linear system of any size. This can be done through Moore-
// Penrose pseudoinverse computation using Singular value decomposition.
// All the other methods provided are for convenience or testing.
//
// An instance of BigMatrix class can grow and shrink, but, when shrunk,
// the matrix only pretends it's smaller (it returns the smaller height
// and width, but the memory usage remains the same). When resizing back
// and forth, no reallocation (and therefore no loosing of cell contents)
// appears. Otherwise (when the matrix needs to grow beyond its current
// cap height or width), the whole matrix is freed and reallocated. The 
// new contents is undefined.
//
// Also, the big matrix uses the simple memory manager: whenever it allocates
// a memory, the memory never goes back tot he OS; it remains 
// in the small pool that is re-used next time the matrix of the same size is 
// needed. Thus, NOTE:
//
// this matrix implementation is optimized for often use throughout the
// program; it behaves badly when the program allocates a lot of matrices
// at once for one-time computation, or when each time the matrix is
// created, it has a new size. It behaves well if the matrices are allocated
// often, are of roughly same size and there are never too many matrices
// needed at once.
//
// the matrix is kept in row-major, so when indexing, the first index
// is the row, the second index is the column
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_BIGMATRIX_H__6148BEAD_11B0_4A41_87F2_61F89675B506__INCLUDED_)
#define AFX_BIGMATRIX_H__6148BEAD_11B0_4A41_87F2_61F89675B506__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

class CBigVector;

/////////////////////////////////////////////////////////////
// class CBigMatrix
//
// The big matrix, a resizable rectangular matrix.
// see also description above.
/////////////////////////////////////////////////////////////
class CBigMatrix  
{
public:
	// creates an empty matrix; needs to be resized
	CBigMatrix();
	
	// copies the size and contents of the matrix
	// the new matrix is allocated with the same capacity as the source
	CBigMatrix(const CBigMatrix& mxThat);

	// prepares a matrix of the specified dimensions
	// the content is undefined
	CBigMatrix(int nHeight, int nWidth);

	~CBigMatrix();

	// copies contents of the matrix
	CBigMatrix& operator = (const CBigMatrix& mxThat);

	// resizes the matrix; doesn't preserve its content, just zeros it out
	void Resize(int nHeight, int nWidth);

	// returns the pointer to a specified row in the matrix.
	// the row is an array of doubles and can be indexed appropriately, e.g.
	// bmMyMatrix[nRow][nColumn]
	double* operator [] (int nRow);
	const double* operator [] (int nRow)const;

	// return dimensions of the matrix
	int GetHeight() const {return m_nSizeHeight;}
	int GetWidth() const {return m_nSizeWidth;}
	
	// reutrn capacity of the matrix, that is, the maximum dimensions to
	// which it can be resized without reallocation
	int GetCapHeight() const {return m_nCapHeight;}
	int GetCapWidth() const {return m_nCapWidth;}


	// composes pseudoinverse of this matrix = (A'*A)^-1*A', where ' is transpose and ^-1 is inverse
	// if there is an entity to divide to that is smaller than the given threshold,
	// the process is aborted and an empty matrix is returned
	// RETURNS: 0 if failed; 1 if the found matrix is VERY stable
	//          something in - between otherwise (lim 0 is unstable, 0.1..1 is stable area)
	double ComposePseudoInverse (CBigMatrix& mxPseudoInverse, double fThreshold)const;


	// composes inverse of this matrix; if in the process we have to divide by
	// fThreshold or less then we return an empty matrix
	void ComposeInverse (CBigMatrix& mxInverse, double fThreshold);

	// sets this matrix to identity
	void SetIdentity();

	// sets this matrix to production of the given arguments (A*B)
	void SetMultiply (CBigMatrix& mxA, CBigMatrix& mxB);

	// sets the matrix to subtraction A-B
	void SetSubtract (CBigMatrix& mxA, CBigMatrix& mxB);

	// returns the square root of the sum of squares of this matrix
	double GetNorm ()const;
	
	// returns the square root of the sum of squares of (this matrix minus identity matrix)
	double GetIdentityError()const;

	// assigns transpose of this matrix to the target matrix
	void Transpose (CBigMatrix& mxT);

	// zeroes this matrix, including the spare cells (if capacity > dimensions)
	void SetZero();

	// orthogonalizes and unifies the matrix
	//void OrthoUnify
protected:
	// Singular Value Decomposition, returns true if succeeded
	// NOTE: converts this == SVD(this) * DIAG(w) * V
	bool SVDecomposition (CBigVector& w, CBigMatrix& V, int nMaxIterations = 30);
	bool SVDecompositionV2 (CBigVector& w, CBigMatrix& V, int nMaxIterations = 30);
	bool SVDecompositionV3 (CBigVector& w, CBigMatrix& V, int nMaxIterations = 30);

	// composes inverse of this matrix; if in the process we have to divide by
	// fThreshold or less then we return an empty matrix
	// SIDE EFFECT: sets this matrix to identity
	void internalComposeInverse (CBigMatrix& mxInverse, double fThreshold);

	// the matrix current capacity
	int m_nCapHeight, m_nCapWidth;
	// the matrix current size; cannot be larger than the capacity in any dimension	
	int m_nSizeHeight, m_nSizeWidth;
	// the matrix cells, row-major order
	double* m_fCells;
};

// operators for matrix-vector multiplication, follow the normal linear algebra rules
// I wouldn't recommend to use this because of it being extremely slow
// provided for testing or usage in non-performance-critical parts of the code only
extern CBigVector operator * (const CBigVector&, const CBigMatrix&);
extern CBigVector operator * (const CBigMatrix&, const CBigVector&);

#endif // !defined(AFX_BIGMATRIX_H__6148BEAD_11B0_4A41_87F2_61F89675B506__INCLUDED_)

