tcovar.cpp

这是个人脸识别程序

			jMax = nelems;
		}
	else
		jMax = nelems;
	for (int j = i+1; j < jMax; j++)
		Dist += pA1[j] * xi * xTemp[j];			// pA1[j] is A(i,j)
	}
Dist *= 2;										// account for lower left triangle elements (they were not included in the above loops)

for (i = 0; i < nelems; i++)					// sum diag elements
	{
	double *const pA1 = pA + (i * iTda) + i;	// same as A(i,i)
	Dist += *pA1 * xTemp[i] * xTemp[i];
	}
return Dist;
}

//-----------------------------------------------------------------------------
// Convert the given matrix, which must be a symmetric matrix, to sparse format.
// We only store entries for the upper triangle, that's why the matrix should be symmetric.
//
// The sparse format is:
//
//		iEntry  Col0  Col1  Col2
//		0       nrows ncols MagicNbr
//		1..n    iRow  iCol  value
//
// The diagonal entries are first, then the rest of the entries.
// All diagonal entries are always stored even if they are 0.
//
// When building the sparse matrix, the lower left triangle elements of A
// are ignored because it's assumed A is symmetrical

void ConvertSymmetricMatToSparseFormat (Mat &A, 			// io
										double MinVal)		// in: clear any element smaller than this
{
int nrows = A.nrows(), ncols = A.ncols();
Mat B(nrows * ncols, 3);		// working array

B(0,0) = nrows;
B(0,1) = ncols;
B(0,2) = SPARSE_SYMMETRIC;

int nEntries = 1;
ASSERT(nrows == ncols);

// Diag entries.
// We always include all diag entries, even if they are 0.
// CopyMatToSparseMat assumes that all diag entries are present.

for (int iDiag = 0; iDiag < nrows; iDiag++)
	{
	B(nEntries, 0) = iDiag;							
	B(nEntries, 1) = iDiag;
	B(nEntries, 2) = A(iDiag, iDiag);
	nEntries++;
	}
// upper triangle entries

for (int iRow = 0; iRow < nrows; iRow++)
	for (int iCol = iRow+1; iCol < ncols; iCol++)	// start at iRow+1 to get upper triangle only
		if (fabs(A(iRow, iCol)) > MinVal)
			{
			B(nEntries, 0) = iRow;
			B(nEntries, 1) = iCol;
			B(nEntries, 2) = A(iRow, iCol);
			nEntries++;
			}

// copy working array B back into A

A.dim(nEntries, 3);
for (int iEntry = 0; iEntry < nEntries; iEntry++)
	{
	A(iEntry, 0) = B(iEntry, 0);
	A(iEntry, 1) = B(iEntry, 1);
	A(iEntry, 2) = B(iEntry, 2);
	}
}

//-----------------------------------------------------------------------------
// This makes assumptions about ordering of elements in A: first the diag elems
// then the upper triangular elems in col major order

static void PrintSparseMat (const Mat &A)
{
int nrows = A(0,0);
int ncols = A(0,1);
int iEntry = 1;
int iEntry1 = nrows+1;

lprintf("\n");

for (int iRow = 0; iRow < nrows; iRow++)
	{
	for (int iCol = 0; iCol < ncols; iCol++)
		if (iEntry < A.nrows() && A(iEntry,0) == iRow && A(iEntry,1) == iCol)			// diag entry?
			{
			lprintf("%g ", A(iEntry, 2));
			iEntry++;
			}
		else if (iEntry1 < A.nrows() && A(iEntry1,0) == iRow && A(iEntry1,1) == iCol)	// non diag entry?
			{
			lprintf("%g ", A(iEntry1, 2));
			iEntry1++;
			}
		else
			lprintf("0 ");

	lprintf("\n");
	}
}

//-----------------------------------------------------------------------------
// this copies the Mat M to the SparseMat S, after resizing S appropriately

void CopyMatToSparseMat (SparseMat &S,  // out
                         const Mat &M)  // in
{
const int nSparseRows = M.nrows();
S.resize(nSparseRows);

for (int i = 0; i < nSparseRows; i++)
    {
    S[i].iRow = short(M(i, 0));     // type convert: double to short
    S[i].iCol = short(M(i, 1));     // ditto
    S[i].Val  = M(i, 2);
    }
}

//-----------------------------------------------------------------------------
static void PrintSparseMat (const SparseMat &A)
{
int nrows = A[0].iRow;
int ncols = A[0].iCol;
int iEntry = 1;
int iEntry1 = nrows+1;

lprintf("\n");

for (int iRow = 0; iRow < nrows; iRow++)
	{
	for (int iCol = 0; iCol < ncols; iCol++)
		{
		if (iEntry < A.size() && A[iEntry].iRow == iRow && A[iEntry].iCol == iCol)			// diag entry?
			{
			lprintf("%g ", A[iEntry].Val);
			iEntry++;
			}
		else if (iEntry1 < A.size() && A[iEntry1].iRow == iRow && A[iEntry1].iCol == iCol)	// non diag entry?
			{
			lprintf("%g ", A[iEntry1].Val);
			iEntry1++;
			}
		else
			lprintf("0 ");
		}
	lprintf("\n");
	}
}

//-----------------------------------------------------------------------------
// Returns true if mat is sparse.
// Can return false positives if A[0,2] happens to equal the constant SPARSE_SYMMETRIC
// But will never return a false positive if ncols !=3

bool fSparseMat (const Mat &A)
{
return A.ncols() == 3 && A(0,2) == SPARSE_SYMMETRIC;
}

//-----------------------------------------------------------------------------
// Return x' * A * x where A is a sparse square symmetric matrix and x is a vec.
//
// This function is equivalent to OneElemToDouble(x.t() * A * x), but is
// optimized for speed and is much faster, if the matrix is sparse enough.
// Efficiency is paramount in this routine.

double Sparse_xAx (const Vec &x, const SparseMat &A)    // in: all args
{
const int nx = x.nelems();
ASSERT(x.ncols() == 1 || x.nrows() == 1);
ASSERT(A[0].iRow == nx);                // first elem is nrows,ncols,SPARSE_SYMMETRIC
ASSERT(A[0].iCol == nx);
ASSERT(A[0].Val == SPARSE_SYMMETRIC);
ASSERT(A.size() > unsigned(nx));        // > not >= because first elem is not mat data

double DiagResult = 0, Result = 0;

// sum diag elements

const double * const px = x.m->data;    // for efficiency, access mat buf directly
int i;
for (i = 0; i < nx; i++)
    {
    DASSERT(A[i+1].Val >= 0);
    DASSERT(A[i+1].iRow == A[i+1].iCol);
    DiagResult += A[i+1].Val * px[i] * px[i];
    }
// sum upper right triangle elements

const int nSparseRows = A.size();
for (i++; i < nSparseRows; i++)
    {
    DASSERT(A[i].iRow != A[i].iCol);
    Result += px[A[i].iRow] * px[A[i].iCol] * A[i].Val;
    }
Result *= 2;                            // incorporate lower left triangle elements

return DiagResult + Result;
}