Subversion Repositories or1k

/***********************************************************************
 
sc.c - creates classifiers from feature vectors of examples, as well as
   classifying example feature vectors.
 
Copyright (C) 1991 Dean Rubine
 
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 1, 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.
 
You should have received a copy of the GNU General Public License
along with this program (in ../COPYING); if not, write to the Free
Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
***********************************************************************/
 
 
#include "bitvector.h"
#include "matrix.h"
#include "util.h"
#include "sc.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "zdebug.h"
 
#define	EPS	(1.0e-6)	/* for singular matrix check */
sClassifier
sNewClassifier()
{
	register sClassifier sc = allocate(1, struct sclassifier);
	sc->nfeatures = -1;
	sc->nclasses = 0;
	sc->classdope = allocate(MAXSCLASSES, sClassDope);
	sc->w = NULL;
	return sc;
}
 
void
sFreeClassifier(sc)
register sClassifier sc;
{
	register int i;
	register sClassDope scd;
 
	for(i = 0; i < sc->nclasses; i++) {
		scd = sc->classdope[i];
		if(scd->name) free(scd->name);
		if(scd->sumcov) FreeMatrix(scd->sumcov);
		if(scd->average) FreeVector(scd->average);
		free(scd);
		if(sc->w && sc->w[i]) FreeVector(sc->w[i]);
	}
	free(sc->classdope);
	if(sc->w) free(sc->w);
	if(sc->cnst) FreeVector(sc->cnst);
	if(sc->invavgcov) FreeMatrix(sc->invavgcov);
	free(sc);
}
 
sClassDope
sClassNameLookup(sc, classname)
register sClassifier sc;
register char *classname;
{
	register int i;
	register sClassDope scd;
	static sClassifier lastsc = NULL;
	static sClassDope lastscd = NULL;
 
	/* quick check for last class name */
	if(lastsc == sc && lastscd != NULL && STREQ(lastscd->name, classname))
		return lastscd;
 
	/* linear search through all classes for name */
	for(i = 0; i < sc->nclasses; i++) {
		scd = sc->classdope[i];
		if(STREQ(scd->name, classname))
			return lastsc = sc, lastscd = scd;
	}
	lastsc = NULL;
	lastscd = NULL;
	return NULL;
}
 
static sClassDope
sAddClass(sc, classname)
register sClassifier sc;
char *classname;
{
	register sClassDope scd;
 
	sc->classdope[sc->nclasses] = scd = allocate(1, struct sclassdope);
	scd->name = scopy(classname);
	scd->number = sc->nclasses;
	scd->nexamples = 0;
	scd->sumcov = NULL;
	++sc->nclasses;
	return scd;
}
 
void
sAddExample(sc, classname, y)
register sClassifier sc;
char *classname;
Vector y;
{
	register sClassDope scd;
	register int i, j;
	double nfv[50];
	double nm1on, recipn;
 
	scd = sClassNameLookup(sc, classname);
	if(scd == NULL) {
/* fprintf(stderr, "sAddExample: calling sAddClass on %s.\n", classname); */
		scd = sAddClass(sc, classname);
	      }
 
	if(sc->nfeatures == -1) {
		sc->nfeatures = NROWS(y);
/*		fprintf(stderr, "sAddExample: setting sc->nfeatures to NROWS(y).\n"); */
	      }
 
	if(scd->nexamples == 0) {
/* 		fprintf(stderr, "sAddExample: allocating  & zeroing scd->average & scd->sumcov.\n"); */
		scd->average = NewVector(sc->nfeatures);
		ZeroVector(scd->average);
		scd->sumcov = NewMatrix(sc->nfeatures, sc->nfeatures);
		ZeroMatrix(scd->sumcov);
 
	}
 
	if(sc->nfeatures != NROWS(y)) {
		PrintVector(y, "sAddExample: funny feature vector nrows!=%d", 
			sc->nfeatures);
		return;
	}
 
	scd->nexamples++;
	nm1on = ((double) scd->nexamples-1)/scd->nexamples;
	recipn = 1.0/scd->nexamples;
 
	/* incrementally update covariance matrix */
        for(i = 0; i < sc->nfeatures; i++)
		nfv[i] = y[i] - scd->average[i];
 
	/* only upper triangular part computed */
        for(i = 0; i < sc->nfeatures; i++)
	   for(j = i; j < sc->nfeatures; j++)
		scd->sumcov[i][j] += nm1on * nfv[i] * nfv[j];
 
	/* incrementally update mean vector */
	for(i = 0; i < sc->nfeatures; i++)
		scd->average[i] = nm1on * scd->average[i] + recipn * y[i];
 
}
 
void
sDoneAdding(sc)
register sClassifier sc;
{
	register int i, j;
	int c;
	int ne, denom;
	double oneoverdenom;
	register Matrix s;
	register Matrix avgcov;
	double det;
	register sClassDope scd;
 
	if(sc->nclasses == 0) {
		error("No classes for adding to classifier");
		return;
	    }
 
	/* Given covariance matrices for each class (* number of examples - 1)
	    compute the average (common) covariance matrix */
 
	avgcov = NewMatrix(sc->nfeatures, sc->nfeatures);
	ZeroMatrix(avgcov);
	ne = 0;
	for(c = 0; c < sc->nclasses; c++) {
		scd = sc->classdope[c];
		ne += scd->nexamples;
		s = scd->sumcov;
		for(i = 0; i < sc->nfeatures; i++)
			for(j = i; j < sc->nfeatures; j++)
				avgcov[i][j] += s[i][j]; 
	}
 
	denom = ne - sc->nclasses;
	if(denom <= 0) {
	    error("Number of classes must be less than number of examples");
	    return;
	}
 
	oneoverdenom = 1.0 / denom;
	for(i = 0; i < sc->nfeatures; i++)
		for(j = i; j < sc->nfeatures; j++)
			avgcov[j][i] = avgcov[i][j] *= oneoverdenom;
 
	Z('a') PrintMatrix(avgcov, "Average Covariance Matrix\n");
	/* invert the avg covariance matrix */
 
	sc->invavgcov = NewMatrix(sc->nfeatures, sc->nfeatures);
	det = InvertMatrix(avgcov, sc->invavgcov);
	if(fabs(det) <= EPS)
		FixClassifier(sc, avgcov);
 
	/* now compute discrimination functions */
	sc->w = allocate(sc->nclasses, Vector);
	sc->cnst = NewVector(sc->nclasses);
	for(c = 0; c < sc->nclasses; c++) {
		scd = sc->classdope[c];
		sc->w[c] = NewVector(sc->nfeatures);
		VectorTimesMatrix(scd->average, sc->invavgcov, sc->w[c]);
		sc->cnst[c] = -0.5 * InnerProduct(sc->w[c], scd->average);
		/* could add log(priorprob class c) to cnst[c] */
	}
 
	FreeMatrix(avgcov);
	return;
}
 
sClassDope
sClassify(sc, fv) {
	return sClassifyAD(sc, fv, NULL, NULL);
}
 
sClassDope
sClassifyAD(sc, fv, ap, dp)
sClassifier sc;
Vector fv;
double *ap;
double *dp;
{
	double disc[MAXSCLASSES];
	register int i, maxclass;
	double denom, exp();
	register sClassDope scd;
	double d;
 
	if(sc->w == NULL) {
		error("%x not a trained classifier", sc);
		return(NULL);
       }
 
	for(i = 0; i < sc->nclasses; i++) {
/* ari */
	  double IP;
	  IP = InnerProduct(sc->w[i], fv);
/*	  fprintf(stderr, "sClassifyAD:  InnerProduct for class %s is %f.\n", sc->classdope[i]->name, IP); */
/*	  fprintf(stderr, "sClassifyAD:  sc->cnst[i] = %f.\n", sc->cnst[i]); */
	  disc[i] = IP + sc->cnst[i];
/*	  fprintf(stderr, "sClassifyAD:  Set disc = %f for class %s.\n", disc[i],sc->classdope[i]->name); */
 
/*	  disc[i] = InnerProduct(sc->w[i], fv) + sc->cnst[i]; */
	}
 
	maxclass = 0;
	for(i = 1; i < sc->nclasses; i++)
		if(disc[i] > disc[maxclass])
			maxclass = i;
 
/* ari */
/* PF_INIT_COS	0	 initial angle (cos)                         */
/* PF_INIT_SIN	1	 initial angle (sin)                         */
/* PF_BB_LEN	2	 length of bounding box diagonal             */
/* PF_BB_TH	3	 angle of bounding box diagonal              */
/* PF_SE_LEN	4	 length between start and end points         */
/* PF_SE_COS	5	 cos of angle between start and end points   */
/* PF_SE_SIN	6	 sin of angle between start and end points   */
/* PF_LEN	7	 arc length of path                          */
/* PF_TH	8	 total angle traversed                       */
/* PF_ATH	9	 sum of abs vals of angles traversed         */
/* PF_SQTH	10	 sum of squares of angles traversed          */
/* PF_DUR	11	 duration of path                            */ 
/* ifndef USE_TIME                                                   */
/* 	NFEATURES	12                                           */
/* else                                                              */
/* 	PF_MAXV		12	   maximum speed                     */
/* 	NFEATURES	13                                           */
/* endif                                                             */
 
/*
* fprintf(stderr, "\nFeature vector:\n");
* fprintf(stderr, "    start cosine      %8.4f    path length       %8.4f\n",
* 	fv[PF_INIT_COS], fv[PF_LEN]);
* fprintf(stderr, "    start sine        %8.4f    total angle       %8.4f\n",
* 	fv[PF_INIT_SIN], fv[PF_TH]);
* fprintf(stderr, "    b.b. length       %8.4f    total abs. angle  %8.4f\n",
* 	fv[PF_BB_LEN], fv[PF_ATH]);
* fprintf(stderr, "    b.b. angle        %8.4f    total sq. angle   %8.4f\n",
* 	fv[PF_BB_TH], fv[PF_SQTH]);
* fprintf(stderr, "    st-end length     %8.4f    duration          %8.4f\n",
* 	fv[PF_SE_LEN], fv[PF_DUR]);
* fprintf(stderr, "    st-end cos        %8.4f\n", fv[PF_SE_COS]);
* fprintf(stderr, "    st-end sin        %8.4f\n", fv[PF_SE_SIN]);
*/
 	ZZ('C') {
		scd = sc->classdope[maxclass];
		PrintVector(fv, "%10.10s  ", scd->name);
		ZZZ('C') {
			for(i = 0; i < sc->nclasses; i++) {
				scd = sc->classdope[i];
				PrintVector(scd->average, "%5.5s %5g ",
					scd->name, disc[i]);
			}
		}
	}
 
	scd = sc->classdope[maxclass];
/* ari */
/* fprintf(stderr,"%s", scd->name); */
/*
   fprintf(stderr,"Stroke identified as %s [ ", scd->name);
   for (i = 0; i < sc->nclasses; i++) {
      if ( (disc[maxclass] - disc[i] < 5.0) && (i != maxclass) )
         fprintf(stderr,"%s ", sc->classdope[i]->name);
   }
   fprintf(stderr,"], ");
*/
	if(ap) {	/* calculate probability of non-ambiguity */
		for(denom = 0, i = 0; i < sc->nclasses; i++)
			/* quick check to avoid computing negligible term */
			if((d = disc[i] - disc[maxclass]) > -7.0)
				denom += exp(d);
		*ap = 1.0 / denom;
	}
 
	if(dp) 	/* calculate distance to mean of chosen class */
		*dp = MahalanobisDistance(fv, scd->average, sc->invavgcov);
 
	return scd;
}
 
/*
 Compute (v-u)' sigma (v-u)
 */
 
double
MahalanobisDistance(v, u, sigma)
register Vector v, u;
register Matrix sigma;
{
	register int i;
	static Vector space;
	double result;
 
	if(space == NULL || NROWS(space) != NROWS(v)) {
		if(space) FreeVector(space);
		space = NewVector(NROWS(v));
	}
	for(i = 0; i < NROWS(v); i++)
		space[i] = v[i] - u[i];
	result =  QuadraticForm(space, sigma);
	return result;
}
 
void
FixClassifier(sc, avgcov)
register sClassifier sc;
Matrix avgcov;
{
	int i;
	double det;
	BitVector bv;
	Matrix m, r;
 
	/* just add the features one by one, discarding any that cause
	   the matrix to be non-invertible */
 
	CLEAR_BIT_VECTOR(bv);
	for(i = 0; i < sc->nfeatures; i++) {
		BIT_SET(i, bv);
		m = SliceMatrix(avgcov, bv, bv);
		r = NewMatrix(NROWS(m), NCOLS(m));
		det = InvertMatrix(m, r);
		if(fabs(det) <= EPS)
			BIT_CLEAR(i, bv);
		FreeMatrix(m);
		FreeMatrix(r);
	}
 
	m = SliceMatrix(avgcov, bv, bv);
	r = NewMatrix(NROWS(m), NCOLS(m));
	det = InvertMatrix(m, r);
	if(fabs(det) <= EPS) {
		error("Can't fix classifier!");
		return;
	    }
	DeSliceMatrix(r, 0.0, bv, bv, sc->invavgcov);
 
	FreeMatrix(m);
	FreeMatrix(r);
 
}
 
void
sDumpClassifier(sc)
register sClassifier sc;
{
	register sClassIndex c;
 
	printf("\n----Classifier %x, %d features:-----\n", (int)sc, sc->nfeatures);
	printf("%d classes: ", sc->nclasses);
	for(c = 0; c < sc->nclasses; c++)
		printf("%s  ", sc->classdope[c]->name);
	printf("Discrimination functions:\n");
	for(c = 0; c < sc->nclasses; c++) {
		PrintVector(sc->w[c], "%s: %g + ", sc->classdope[c]->name, sc->cnst[c]);
		printf("Mean vectors:\n");
		PrintVector(sc->classdope[c]->average, "%s: ", sc->classdope[c]->name);
	    }
	if( sc->invavgcov != NULL ) {
	    PrintMatrix(sc->invavgcov, "Inverse average covariance matrix:\n");
	}
	printf("\n---------\n\n");
}
 
void
sWrite(outfile, sc)
FILE *outfile;
sClassifier sc;
{
	int i;
	register sClassDope scd;
 
	fprintf(outfile, "%d classes\n", sc->nclasses);
	for(i = 0; i < sc->nclasses; i++) {
		scd = sc->classdope[i];
		fprintf(outfile, "%s\n", scd->name);
	}
	for(i = 0; i < sc->nclasses; i++) {
		scd = sc->classdope[i];
		OutputVector(outfile, scd->average);
		OutputMatrix(outfile, scd->sumcov);
		OutputVector(outfile, sc->w[i]);
	}
	OutputVector(outfile, sc->cnst);
	OutputMatrix(outfile, sc->invavgcov);
}
 
sClassifier
sRead(infile)
FILE *infile;
{
	int i, n;
	register sClassifier sc;
	register sClassDope scd;
	char buf[100];
 
 
	Z('a') printf("Reading classifier \n");
	sc = sNewClassifier();
	fgets(buf, 100, infile);
	if(sscanf(buf, "%d", &n) != 1) {
	    error("Input error in classifier file");
	    sFreeClassifier(sc);
	    return(NULL);
	}
	Z('a') printf("  %d classes \n", n);
	for(i = 0; i < n; i++) {
		fscanf(infile, "%s", buf);
		scd = sAddClass(sc, buf);
		Z('a') printf("  %s \n", scd->name);
	}
	sc->w = allocate(sc->nclasses, Vector);
	for(i = 0; i < sc->nclasses; i++) {
		scd = sc->classdope[i];
		scd->average = InputVector(infile);
		scd->sumcov = InputMatrix(infile);
		sc->w[i] = InputVector(infile);
	}
	sc->cnst = InputVector(infile);
	sc->invavgcov = InputMatrix(infile);
	Z('a') printf("\n");
	return sc;
}
 
 
void
sDistances(sc, nclosest)
register sClassifier sc;
{
	register Matrix d = NewMatrix(sc->nclasses, sc->nclasses);
	register int i, j;
	double min, max = 0;
	int n, mi, mj;
 
	printf("-----------\n");
	printf("Computing %d closest pairs of classes\n", nclosest);
	for(i = 0; i < NROWS(d); i++) {
		for(j = i+1; j < NCOLS(d); j++) {
			d[i][j] = MahalanobisDistance(
						sc->classdope[i]->average,
						sc->classdope[j]->average,
						sc->invavgcov);
			if(d[i][j] > max) max = d[i][j];
		}
	}
 
	for(n = 1; n <= nclosest; n++) {
		min = max;
		mi = mj = -1;
		for(i = 0; i < NROWS(d); i++) {
			for(j = i+1; j < NCOLS(d); j++) {
				if(d[i][j] < min)
					min = d[mi=i][mj=j];
			}
		}
		if(mi == -1)
			break;
 
		printf("%2d) %10.10s to %10.10s d=%g nstd=%g\n",
			n,
			sc->classdope[mi]->name,
			sc->classdope[mj]->name,
			d[mi][mj],
			sqrt(d[mi][mj]));
 
		d[mi][mj] = max+1;
	}
	printf("-----------\n");
	FreeMatrix(d);
}