#include <string.h>#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 */sClassifiersNewClassifier(){    register sClassifier sc = allocate(1, struct sclassifier);    sc->nfeatures = -1;    sc->nclasses = 0;    sc->classdope = allocate(MAXSCLASSES, sClassDope);    sc->w = NULL;    return sc;}voidsFreeClassifier(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);}sClassDopesClassNameLookup(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 sClassDopesAddClass(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;}voidsAddExample(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];}voidsDoneAdding(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;}sClassDopesClassify(sc, fv){    return sClassifyAD(sc, fv, NULL, NULL);}sClassDopesClassifyAD(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) */doubleMahalanobisDistance(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;}voidFixClassifier(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);}voidsDumpClassifier(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");}voidsWrite(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);}sClassifiersRead(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;}voidsDistances(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);}