PolynomialLinearRegression.cpp

#include "StdAfx.h"
#include "PolynomialLinearRegression.h"

CPolynomialLinearRegression::CPolynomialLinearRegression(LR_CountT nRegressors)
{
	m_pSumXn = NULL;
	ReInit(nRegressors);
}
void CPolynomialLinearRegression::ReInit(LR_CountT nRegressors)
{
	LR_CountT nSize;
	LR_CountT nSumXnCnt;
	LR_CountT nSumXnYCnt;
	LR_CountT nResultsCnt;
	if( nRegressors < 1 )
	{
		nRegressors = 1;
	}
	nSumXnCnt = (nRegressors * 2) + 1;
	nSumXnYCnt = nRegressors + 1;
	nResultsCnt = nSumXnYCnt;
	nSize = nSumXnCnt + nSumXnYCnt + nResultsCnt;
	do
	{
		if( m_pSumXn != NULL )
		{
			if( m_nRegressors == nRegressors )
			{
				break;
			}
			delete []m_pSumXn;
		}
		m_pSumXn = new LR_DataT[nSize];
		m_pSumXnY = &m_pSumXn[nSumXnCnt];
		m_pResults = &m_pSumXnY[nSumXnYCnt];
	}while(0);
	m_nRegressors = nRegressors;
	
	// for( LR_CountT i = 0 ; i < nSize ; i++ )
	// {
		// m_pSumXn[i] = 0.0;
	// }
	memset(m_pSumXn , 0 , nSize*sizeof(LR_DataT));

	m_bIsRefresh = FALSE;
}

CPolynomialLinearRegression::~CPolynomialLinearRegression()
{
	delete []m_pSumXn;
	m_pSumXn = NULL;
}

BOOL CPolynomialLinearRegression::AddXY(LR_DataT x , LR_DataT y)
{
	LR_CountT i;
	LR_CountT nSumXnCntD;
	LR_CountT nSumXnYCnt;
	LR_DataT Xn;
	nSumXnCntD = m_nRegressors * 2;
	nSumXnYCnt = m_nRegressors + 1;
	m_pSumXn[0] += 1.0;
	m_pSumXnY[0] += y;
	Xn = x;
	for( i = 1 ; i < nSumXnCntD ; i++ )
	{
		if( i < nSumXnYCnt )
		{
			m_pSumXnY[i] += Xn * y;
		}
		m_pSumXn[i] += Xn;
		Xn *= x;
	}
	m_pSumXn[i] += Xn;
	
	m_bIsRefresh = FALSE;
	return TRUE;
}

LR_DataT* CPolynomialLinearRegression::GetResults()
{
	if( !m_bIsRefresh )
	{
		if( !Calculate() )
		{
			return NULL;
		}
	}
	return m_pResults;
}

BOOL CPolynomialLinearRegression::IsValid()
{
	return (m_pSumXn[0] >= (LR_DataT)(m_nRegressors + 1));
}

BOOL CPolynomialLinearRegression::Calculate()
{
	LR_CountT i , j , k;
	LR_CountT nRows;
	LR_CountT nColumns;
	LR_CountT nSumXnCnt;
	LR_CountT iMaxIdx;
	LR_DataT* pMatrix;
	LR_DataT** ppIndex;
	LR_DataT* pMaxRow;
	LR_DataT Val;
	
	nRows = m_nRegressors + 1;
	nColumns = m_nRegressors + 2;
	nSumXnCnt = (m_nRegressors * 2) + 1;
	
	pMatrix = new LR_DataT[nColumns * nRows];
	ppIndex = new LR_DataT*[nRows];
	for( i = 0 ; i < nRows ; i++ )
	{
		ppIndex[i] = &pMatrix[i*nColumns];
	}
	m_bIsRefresh = FALSE;

	/*
	E:Sigma
	num   Ex1   Ex2   Ex3   Ex4   Ey
	Ex1   Ex2   Ex3   Ex4   Ex5   Ex1*y
	Ex2   Ex3   Ex4   Ex5   Ex6   Ex2*y
	Ex3   Ex4   Ex5   Ex6   Ex7   Ex3*y
	Ex4   Ex5   Ex6   Ex7   Ex8   Ex4*y
	*/
	// Filling polynomial matrix
	for( i = 0 ; i < nSumXnCnt ; i++ )
	{
		if( i < nRows ) // nSumXnYCnt
		{
			for( j = 0 ; j <= i ; j++ )
			{
				pMatrix[(j*nColumns) + (i-j)] = m_pSumXn[i]; // Row = j , Col = i - j
			}
			pMatrix[(i*nColumns) + (nColumns-1)] = m_pSumXnY[i]; // Row = i , Col = nColumns - 1
		}
		else
		{
			for( j = (i-nRows)+1 ; j < nRows ; j++ )
			{
				pMatrix[(j*nColumns) + (i-j)] = m_pSumXn[i]; // Row = j , Col = i - j
			}
		}
	}

	//// Output
	//for( i = 0 ; i < nRows ; i++ )
	//{
	//	for( k = 0 ; k < nColumns ; k++ )
	//	{
	//		CString str;
	//		str.Format(L"%d,%d : %g" , i , k , (double)pMatrix[i*nColumns+k]);
	//		TRACE(str);
	//	}
	//}

	// Gauss elimination
	for( i = 0 ; i < nRows ; i++ )
	{
		Val = LR_ABS(ppIndex[i][i]);
		iMaxIdx = i;
		for( j = i+1 ; j < nRows ; j++ )
		{
			if( Val < LR_ABS(ppIndex[j][i]) )
			{
				Val = LR_ABS(ppIndex[j][i]);
				iMaxIdx = j;
			}
		}
		if( iMaxIdx != i )
		{
			pMaxRow = ppIndex[i];
			ppIndex[i] = ppIndex[iMaxIdx];
			ppIndex[iMaxIdx] = pMaxRow;
		}
		if( Val <= LR_MinValue )
		{
			goto PLR_Calculate_ErrorEnd;
		}
		for( j = i+1 ; j < nRows ; j++ )
		{
			Val = ppIndex[j][i] / ppIndex[i][i];
			for( k = i+1 ; k <= nRows ; k++ )
			{
				ppIndex[j][k] -= ppIndex[i][k] * Val;
			}
			ppIndex[j][i] = 0.0;
		}
	}

	// Polynomial solving
	for( i = nRows ; i != 0 ; )
	{
		i--;
		m_pResults[i] = ppIndex[i][nRows] / ppIndex[i][i]; // nRows : nColumns - 1 (Y)
		if( LR_ABS(ppIndex[i][i]) <= LR_MinValue ) // Save extreme values : m_pResults[i]
		{
			goto PLR_Calculate_ErrorEnd;
		}
		for( j = i ; j != 0 ; )
		{
			j--;
			ppIndex[j][nRows] -= ppIndex[j][i] * m_pResults[i]; // nRows : nColumns - 1 (Y)
			ppIndex[j][i] = 0.0;
		}
	}
	m_bIsRefresh = TRUE;
PLR_Calculate_ErrorEnd:
	delete []pMatrix;
	delete ppIndex;
	return m_bIsRefresh;
}