tcovar.cpp

这是个人脸识别程序

// $masm\tcovar.cpp 1.5 milbo$ routines for dealing with "trimmed" covariance matrices
// Warning: this is raw research code -- expect it to be quite messy.
//
// Some notes on the parameter nTrimCovar
//
// nTrimCovar is used during training and only for 2d profiles.
//
// If nonzero, trim the covariance matrix so that each profile element can be influenced
// only by other profile elements that are no more than nTrimCovar-1 pixels away.
// Make nTrimCovar=1 to use only only the diagonal elems (i.e variances, no covariances).
// Make nTrimCovar=0 to use the entire covariance matrix as is (no trimmed distance processing)
// Max nTrimCovar is 1+(nProfWidth2d-1)/2
//
// The main reason we use nTrimCovar is for efficiency -- it reduces the number
// of operations needed to calculate xAx.  The current implementation speeds BioId searches
// up by about 100ms when nTrimCovar==5 and nProfWidth2d==11.
//
// milbo jul06 Gardens
//-----------------------------------------------------------------------------
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// A copy of the GNU General Public License is available at
// http://www.r-project.org/Licenses/
//-----------------------------------------------------------------------------

#include "all.hpp"

//-----------------------------------------------------------------------------
// Trim the covariance matrix Covar so that each profile element can be influenced
// only by other profile elements that are no more than nDist pixels away
//
// If you want to use all elements then set nDist=nProfWidth2d and
// Covar wil be unchanged.  Although then it is a waste of time to call this routine.
//
// Example: if nProfWidth2d=11 and nDist=5, then the first part
// the covar matrix becomes (where u means used, . means unused i.e. set to 0):
// Notes: the lower triangle of blocks isn't filled in
//        there are 6=5+1 in total blocks used for each level
//
//             Block0     Block1     Block2     Block3     Block4     Block5     Block6     Blocks7-10...
//  Block0  .: uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.......................... etc.
//          1: uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu......................... etc.
//          2: uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu........................ etc.
//          3: uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu....................... etc.
//          4: uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu...................... etc.
//          5: uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu..................... etc.
//          6: .uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu..................... etc.
//          7: ..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..................... etc.
//          8: ...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu..................... etc.
//          9: ....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu..................... etc.
//         10: .....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu..................... etc.
//  Block1 11: ...........uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu............... etc.
//         12: ...........uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu.............. etc.
//         13: ...........uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu............. etc.
//         14: ...........uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu............ etc.
//         15: ...........uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu........... etc.
//         16: ...........uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu.......... etc.
//         17: ............uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.uuuuuuuuuu.......... etc.
//         18: .............uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu.......... etc.
//         19: ..............uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu.......... etc.
//         20: ...............uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu.......... etc.
//         21: ................uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.......... etc.
//  Block2 22: ......................uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.....uuuuuu.... etc.
//         23: ......................uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu....uuuuuuu... etc.
//         24: ......................uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu...uuuuuuuu.. etc.
//         25: ......................uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu..uuuuuuuuu. etc.
// continues...

static void TrimCovar (Mat &Covar, int nDist, unsigned ProfSpec)
{
DASSERT(ProfSpec & PROF_2d);
DASSERT(nDist > 0);

int width = nGetProfWidth(Covar.ncols(), ProfSpec);
const int nDist1 = nDist - 1;		// so parameter nDist=1 corresponds only to diag elems
DASSERT(nDist1 <= (width-1)/ 2);

Mat NewCovar(Covar.nrows(), Covar.ncols());
NewCovar.fill(0);

// upper triangle blocks (each block is width x width)

for (int iyBlock = 0; iyBlock < width; iyBlock++)
	for (int ixBlock = iyBlock+1; ixBlock < width; ixBlock++)
		if (abs(ixBlock - iyBlock) <= nDist1)
			{
			const int iyMin = iyBlock * width;
			const int iyMax = iyMin + width;
			for (int iy = iyMin; iy < iyMax; iy++)
				{
				const int ixMin = __max(0, iy % width - nDist1) + ixBlock * width;
				const int ixMax = __min(width-1, iy % width + nDist1) + ixBlock * width;
				for (int ix = ixMin; ix <= ixMax; ix++)
					NewCovar(iy, ix) = Covar(iy, ix);
				}
			}

// diagonal blocks (done separately for efficiency -- reduces the number of compares in inner loops)

for (int iBlock = 0; iBlock < width; iBlock++)
	{
	const int iyMin = iBlock * width;
	const int iyMax = iyMin + width;
	for (int iy = iyMin; iy < iyMax; iy++)
		{
		const int ixMin = __max(0, iy % width - nDist1) + iBlock * width;
		const int ixMax = __min(width-1, iy % width + nDist1) + iBlock * width;
		for (int ix = ixMin; ix <= ixMax; ix++)
			NewCovar(iy, ix) = Covar(iy, ix);
		}
	}
Covar = NewCovar;

for (int i = 0; i < Covar.nrows(); i++)			// make symmetrical again (so we can check it later for pos definiteness)
	for (int j = 0; j < i; j++)
		Covar(i,j) = Covar(j,i);
}

//-----------------------------------------------------------------------------
// This forces A to be positive definite by doing the following:
//   Do spectral decomposition: A = Q D Q'  where D is a diagonal matrix of (sorted) eig values
//   Set small or negative elements of D to a small value
//   Reconstruct A using the new D in A = Q D Q'

static void ForcePosDef (Mat &A, 			// io:
	 				 	 double MaxRatio)	// in:
{
Vec EigVals(A.nrows());
Vec EigVecs = GetEigsForSymMat(A, EigVals);			// returns sorted eigs
Mat Diag(A.nrows(), A.nrows());						// all zeroes to start off with

double Min = A.trace() / (A.nrows() * MaxRatio);
DASSERT(Min > 0);
for (int i = 0; i < EigVals.nelems(); i++)
	if (EigVals(i) > Min)
		Diag(i, i) = EigVals(i);
	else
		Diag(i, i) = Min;

A = EigVecs * Diag * EigVecs.t();
}

//-----------------------------------------------------------------------------
// Trim the matrix and force it to be positive definite
// Returns true if the covar matrix is pos def
// TODO this is not very efficient code, it doesn't really matter

bool fIterativeTrimMat (Mat &Covar, int nTrimCovar, unsigned ProfSpec)
{
int iMaxIters = 0;			// safety counter
bool fPosDef;

TrimCovar(Covar, nTrimCovar, ProfSpec);

while (!(fPosDef = fPosDef1(Covar, 0)) && ++iMaxIters < 2)	// 2 gives as good results as any higher number
	{
	ForcePosDef(Covar, 1e3);								// force matrix positive def
	TrimCovar(Covar, nTrimCovar, ProfSpec);
	}
return fPosDef;
}

//-----------------------------------------------------------------------------
// Return x * A * x where A is a trimmed square symmetric matrix and x is a vector.
//
// This is only needed for old format .asm files
// Efficiency is paramount in this routine.
//
// TODO Not sure if the following comments are still correct:
// We do the processing in two phases, to keep the amount of processing in GetProfDist to a mimimum.
// Phase1: call TrimCovar, once for all covar matrices to be trimmed.
// Phase2: in each call to GetProfDist, we skip all-zero blocks in the covar matrix
// If we didn't do Phase2 we would get the same results -- i.e. Phase1 trims the covar matrix
// but Phase2 reduces the processing time.
// Note that in Phase2 we skip all-zero blocks but don't skip over individual 0 elements
// in partially zero blocks -- to do so would be less time efficient.

double Trimmed_xAx (const Vec &x, const Mat &A, unsigned SubProfSpec, int nTrimCovar)
{
double xTemp[MAX_TEMP_ARRAY_LEN];				// use an array rather than a Vec because access is faster
const int nelems = x.ncols();
for (int i = 0; i < nelems; i++)
	xTemp[i] = x(0,i);							// use two indices to access a vector because it's faster

double *const pA = A.m->data;					// access matrix data directly for speed
const int iTda = A.m->tda;
double Dist = 0;

int nProfWidth;
int jMax1;										// will only be used if nTrimCovar != 0
if (nTrimCovar)
	{
	DASSERT(nTrimCovar < nelems);
	nProfWidth = nGetProfWidth(nelems, SubProfSpec);
	jMax1 = nProfWidth * nTrimCovar;
	}
for (i = 0; i < nelems; i++)					// sum off-diag elements, upper right triangle only
	{
	const double xi = xTemp[i];
	double *pA1 = pA + (i * iTda);
	int jMax;
	if (nTrimCovar && (SubProfSpec & PROF_2d))	// two dimensional profile?
		{										// if so, ignore all-zero elements in extreme upper right triangle
		jMax = i + jMax1;
		if (jMax > nelems)