/*
File name: parametric.c
Created by: Ljubomir Buturovic
Created: 08/03/2004
Purpose: implementation of parametric classifiers (linear and quadratic).
*/
/*
Copyright 2004 Ljubomir J. Buturovic
Permission is hereby granted, free of charge, to any person
obtaining a copy of this software and associated documentation files
(the "Software"), to deal in the Software without restriction,
including without limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of the Software,
and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
static char rcsid[] = "$Id: parametric.c,v 1.11 2005/02/23 05:48:00 ljubomir Exp $";
#include <stdlib.h>
#include <errno.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include "parametric.h"
#include "dataset.h"
#include "lin.h"
#include "lau.h"
#include "lmat.h"
#include "lerr.h"
/*
Return weighted class-conditional estimate of the covariance matrix
for 'dset'. The estimate is the sum of class-conditional covariance
matrices of individual classes, weighted by the class probabilities.
Return NULL in case of error and set errno.
*/
static float **weighted_est(struct dataset *dset)
{
int i;
int j;
int k;
int status;
int d;
int offset;
float factor;
float **sigma;
float **csigma; /* class-conditional covariance matrix */
sigma = (float **) 0;
if (dset)
{
d = dset->d;
sigma = fmx_alloc(d, d);
if (sigma)
{
fmx_set(sigma, d, d, 0.0);
status = 0;
offset = 0;
for (i = 0; (i < dset->c) && !status; i++)
{
csigma = cest(&dset->x[offset], d, dset->nd[i], COVARIANCE_MATRIX);
if (csigma)
{
factor = dset->nd[i];
factor = factor/dset->nv;
for (j = 0; j < dset->d; j++)
for (k = 0; k < dset->d; k++)
sigma[j][k] += factor*csigma[j][k];
mx_free((void **) csigma, d);
offset += dset->nd[i];
}
else
{
status = -1;
sigma = (float **) mx_free((void **) sigma, d);
}
}
}
}
return sigma;
}
/*
Return parametric linear classifier trained on 'dset'. The
classifier is a matrix c by d+1, where 'c' is number of classes and
'd' is dimension of input vectors.
In case of error, return NULL and set errc to the error code. Error
codes are EINVAL for bad input arguments, malloc() error codes and
LERR_SINGCOV for singular covariance matrix.
Uses equations (60-62) from R. O. Duda, P. E. Hart and D. G. Stork,
Pattern Classification, Second Edition, John Wiley & Sons, Inc.,
2001.
*/
float **lin_learn(int mode, struct dataset *dset, int *errc)
{
int i;
int j;
int k;
int d;
int c;
int offset;
float fnv;
float *mean;
float **wmx;
float **sigma; /* estimate of cumulative covariance matrix */
float **sig; /* inverse of sigma */
double s;
wmx = (float **) 0;
if (dset && errc)
{
*errc = 0;
d = dset->d;
c = dset->c;
fnv = dset->nv;
if (mode == CUMULATIVE_COV)
sigma = cest(dset->x, d, dset->nv, COVARIANCE_MATRIX);
else
sigma = weighted_est(dset);
if (sigma)
{
sig = fmx_inv(sigma, d, (float *) 0, (float *) 0, errc);
mx_free((void **) sigma, d);
if (sig)
{
wmx = fmx_alloc(c, d+1);
if (wmx)
{
offset = 0;
for (i = 0; (i < c) && !*errc; i++)
{
mean = mest(&dset->x[offset], d, dset->nd[i]);
if (mean)
{
offset += dset->nd[i];
for (j = 0; j < d; j++)
{
s = 0;
for (k = 0; k < d; k++)
s += sig[j][k]*mean[k];
wmx[i][j] = s;
}
s = 0;
for (j = 0; j < d; j++)
s += mean[j]*wmx[i][j];
wmx[i][d] = -0.5*s+log(dset->nd[i]/fnv);
free(mean);
}
else
{
*errc = errno;
wmx = (float **) mx_free((void **) wmx, c);
}
}
}
else
*errc = errno;
mx_free((void **) sig, d);
}
}
else
*errc = errno;
}
else if (errc)
*errc = EINVAL;
return wmx;
}
void qmodel_free(struct qmodel *model, int c)
{
if (model)
{
mx_free((void **) model->wmx, c);
mx_free((void **) model->sigma, c);
vx_free(model);
}
}
/*
Return partial parametric quadratic classifier trained on
'dset'. The partial classifier is a matrix c by d+1, where 'c' is
number of classes and 'd' is dimension of input vectors.
The function uses equations (66-69) from R. O. Duda, P. E. Hart and
D. G. Stork, Pattern Classification, Second Edition, John Wiley &
Sons, Inc., 2001, Section 2.6.3. According to these equations, the
quadratic classifier consists of a quadratic and linear parts. The
quadratic part is trivial; it is simply -0.5 times the inverted
covariance matrix for the corresponding class. Therefore the
function does not return it, but instead merely computes the
class-conditional inverted covariance matrices and stores them in
dset->sigma. The returned array is the linear part of the
classifier.
In case of error, return NULL and set errc to the error code. Error
codes are EINVAL for bad input arguments, malloc() error codes and
LERR_SINGCOV for singular covariance matrix.
*/
struct qmodel *pqc_learn(struct dataset *dset, int *errc)
{
int i;
int d;
int c;
int offset;
int status;
int length;
float fnv;
float det;
double s;
float *mean;
float **wmx;
struct qmodel *pqc;
wmx = (float **) 0;
c = 0;
length = 0;
pqc = (struct qmodel *) 0;
if (dset && errc)
{
*errc = 0;
d = dset->d;
c = dset->c;
length = d*(d+1)/2;
pqc = malloc(sizeof(struct qmodel));
if (pqc)
{
/*
Compute class-conditional inverse covariance matrices.
*/
status = dataset_sigma(dset, errc);
if (!status)
{
wmx = fmx_alloc(c, d+1);
if (wmx)
{
offset = 0;
for (i = 0; (i < c) && !*errc; i++)
{
fnv = dset->nd[i];
mean = mest(&dset->x[offset], d, dset->nd[i]);
if (mean)
{
fvec_smx(mean, dset->sigma[i], d, wmx[i]);
/*
Bias (constant) term.
*/
s = fvec_dot(wmx[i], mean, d, errc);
det = log(dset->det[2*i])+dset->det[2*i+1];
wmx[i][d] = -0.5*(s+det)+log(dset->nd[i]/fnv);
free(mean);
offset += dset->nd[i];
}
else
{
*errc = errno;
wmx = (float **) mx_free((void **) wmx, c);
free(pqc);
pqc = (struct qmodel *) 0;
}
}
}
else
{
*errc = errno;
free(pqc);
pqc = (struct qmodel *) 0;
}
}
else
{
free(pqc);
pqc = (struct qmodel *) 0;
}
}
else
*errc = errno;
}
else if (errc)
*errc = EINVAL;
if (errc && !*errc)
{
pqc->wmx = wmx;
pqc->sigma = fmx_clone(dset->sigma, c, length);
if (!pqc->sigma)
{
if (errc)
*errc = errno;
mx_free((void **) wmx, c);
free(pqc);
pqc = (struct qmodel *) 0;
}
}
return pqc;
}
/*
Save parametric quadratic classifier defined by 'qcx' in
'fname'. 'd' is dimension of the problem, 'c' is number of classes.
The classifier is saved in the following format: for each class, the
linear portion of the classifier is saved first (which is a vector
of length d+1), followed by the inverse of the class-conditional
covariance matrix. To perform classification, use equations (66-69)
from Duda, Hart & Stork.
Return 0 for success, error code for failure.
*/
int pqc_save(struct qmodel *qcx, int d, int c, char *fname)
{
int i;
int j;
int length;
int errc;
int status;
FILE *fptr;
errc = 0;
status = 0;
if (qcx && fname)
{
length = d*(d+1)/2;
fptr = fopen(fname, "w");
if (fptr)
{
for (i = 0; (i < c) && (status >= 0); i++)
{
for (j = 0; (j <= d) && (status >= 0); j++)
status = fprintf(fptr, "%g\t", qcx->wmx[i][j]);
if (status >= 0)
{
/*
Save the inverse covariance matrix in Symmetric
Storage Mode.
*/
for (j = 0; (j < length) && (status >= 0); j++)
{
if (j == length-1)
status = fprintf(fptr, "%g\n", qcx->sigma[i][j]);
else
status = fprintf(fptr, "%g\t", qcx->sigma[i][j]);
}
}
}
status = fclose(fptr);
if (!status)
errc = errno;
/*
Try to close the file even if there was an error.
*/
if (errc)
fclose(fptr);
}
else
errc = errno;
}
else
errc = EINVAL;
return errc;
}
/*
Append parametric quadratic classifier defined by 'qcx' to 'fptr',
in 'bagging' format. The first column is total number of models
'nmodels, followd by model index 'mdx' and model weight 'weight'.
Return 0 for success, error code for failure.
*/
int pqc_write(FILE *fptr, int nmodels, struct qmodel *qcx, int d, int c, int mdx,
float weight)
{
int i;
int j;
int length;
int errc;
int status;
errc = 0;
status = 0;
if (qcx && fptr)
{
length = d*(d+1)/2;
for (i = 0; (i < c) && (status >= 0); i++)
{
fprintf(fptr, "%d\t%d\t%g\t", nmodels, mdx, weight);
for (j = 0; (j <= d) && (status >= 0); j++)
status = fprintf(fptr, "%g\t", qcx->wmx[i][j]);
if (status >= 0)
{
/*
Save the quadratic term coefficients, which are in
fact inverse class-conditional covariance matrix in
Symmetric Storage Mode.
*/
for (j = 0; (j < length) && (status >= 0); j++)
{
if (j == length-1)
status = fprintf(fptr, "%g\n", qcx->sigma[i][j]);
else
status = fprintf(fptr, "%g\t", qcx->sigma[i][j]);
}
}
}
}
else
errc = EINVAL;
return errc;
}
/*
Memory allocator for pqc_load_models(). Return -1 and set errno in
case of failure.
*/
static int pqc_load_models_alloc(struct qmodel ***models, float **wts, char **line,
int model_cnt, int c, int ad, int length,
int llen)
{
int i;
int status;
struct qmodel **qcx;
status = 0;
qcx = malloc(model_cnt*sizeof(struct qmodel *));
if (qcx)
{
for (i = 0; (i < model_cnt) && !status; i++)
{
qcx[i] = malloc(sizeof(struct qmodel));
if (qcx[i])
{
qcx[i]->wmx = fmx_alloc(c, ad);
if (qcx[i]->wmx)
{
qcx[i]->sigma = fmx_alloc(c, length);
if (!qcx[i]->sigma)
status = -1;
}
else
status = -1;
}
else
status = -1;
}
if (!status)
{
*wts = malloc(model_cnt*sizeof(float));
if (*wts)
{
*line = malloc(llen+2);
if (*line)
{
if (model_cnt == 1)
*wts[0] = 1.0;
}
else
status = -1;
}
else
status = -1;
}
else
status = -1;
}
else
status = -1;
if (!status)
*models = qcx;
return status;
}
/*
Load a committee of quadratic classifiers from 'fname'. The function
returns an array of pointers to 'qmodel' structures, one structure
per model. Each struct defines a parametric quadratic classifier.
The function automatically recognizes two types of model files:
single model (MODEL_SINGLE) or multiple model (MODEL_MULTI). Single
model file contains one linear classifier, while multiple model
contains multiple classifiers.
Multiple model file has total number of models in the first column,
followed by the model ID and model weight. The rest of information
is the same as in single model file.
The function sets the model type file in 'type'. For multiple model
file, the number of models is returnes in 'nmodels', and model
weights are returned in 'weights'.
Each element's qmodel->wmx is a 'rows' by 'columns' matrix.
The quadratic coefficients (qmodel->sigma) are simply the inverse
class-conditional covariance matrices. Note that for the actual
classification, the quadratic part has to be multiplied by -0.5.
In case of failure, return NULL and set 'errc'. Possible errors are
memory allocation and file access errors.
*/
struct qmodel **pqc_load_models(char *fname, int *type, int *nmodels, float **weights,
int *rows, int *columns, int *errc)
{
int status;
int i;
int ntk;
int lc;
int llen;
int nmd;
int lcntr;
int fcol;
int model_cnt;
int length;
int offset;
int c;
int d;
int first_line;
float fad;
char *line;
float *wts;
char **tokens;
struct qmodel **models;
FILE *fptr;
if (fname && errc)
{
models = (struct qmodel **) 0;
*errc = 0;
lc = file_info(fname, &llen, &ntk, '\0');
if ((lc != -1) && (ntk > 0))
{
*type = detect_model(fname, &model_cnt);
if (*type != -1)
{
fptr = fopen(fname, "r");
if (fptr)
{
c = lc/model_cnt;
if (rows)
*rows = c;
if (*type == MODEL_MULTI)
fcol = ntk-3;
else
fcol = ntk;
/*
Magic formula to calculate the number of features
d of the data set which was used to build this
quadratic classifier. See R & D, vol. V, page
10/09/2004. The number of columns in the linear
component of the classifier equals d+1.
*/
fad = 8*fcol+1;
fad = 0.5*(sqrt(fad)-3);
d = fad;
if (columns)
*columns = d+1;
length = (d+1)*d/2;
status = pqc_load_models_alloc(&models, &wts, &line, model_cnt,
c, d+1, length, llen);
if (!status)
{
nmd = 0;
lcntr = 0;
first_line = 1;
while (fgets(line, llen+2, fptr) && !status)
{
tokens = str_tokenize(line, WHITESPACE);
if (tokens)
{
offset = 0;
if (*type == MODEL_MULTI)
{
wts[nmd] = atof(tokens[1]);
offset = 3;
}
for (i = 0; i < d+1; i++)
models[nmd]->wmx[lcntr][i] = atof(tokens[i+offset]);
offset += d+1;
for (i = 0; i < length; i++)
models[nmd]->sigma[lcntr][i] = atof(tokens[i+offset]);
if ((lcntr == c-1) && (nmd < model_cnt-1))
{
/*
Next model.
*/
nmd++;
lcntr = 0;
}
else
lcntr++;
str_free(tokens);
}
else
status = -1;
}
}
else
*errc = errno; /* alloc() failure */
status = fclose(fptr);
if (status == -1)
{
for (i = 0; i < model_cnt; i++)
mx_free((void **) models[i], c);
models = vx_free(models);
*errc = errno;
}
else if (weights)
*weights = wts;
}
else
*errc = errno;
}
else
*errc = errno;
}
else
*errc = errno;
if (nmodels)
*nmodels = model_cnt;
}
else if (errc)
*errc = EINVAL;
return models;
}
/*
Internal function for quadratic classifier. Computes the value of
the parametric quadratic classifier discriminant function.
*/
static float pqc_compute(float *model, float *sigma, int d, float *vector, float *fv)
{
float gx;
gx = fvec_dot(model, vector, d, (int *) 0);
gx += model[d];
/*
Add quadratic term.
*/
fvec_smx(vector, sigma, d, fv);
gx += -0.5*fvec_dot(vector, fv, d, (int *) 0);
return gx;
}
/*
Return predicted classification of 'vector' by parametric quadratic
classifer stored in 'model'. It is assumed that the length of
'vector' is ad-1.
Return -1 in case of bad arguments or malloc() error and set
errno. Otherwise, the result is in the [0..c-1] interval.
*/
int pqc_predict(struct qmodel *model, int c, int ad, float *vector)
{
int i;
int d;
int predicted_class;
float gx;
float gmax;
float *fv;
predicted_class = -1;
gmax = -FLT_MAX;
d = ad-1;
if (model && (c > 0) && (ad > 0) && vector)
{
fv = malloc(d*sizeof(float));
if (fv)
{
for (i = 0; i < c; i++)
{
gx = pqc_compute(model->wmx[i], model->sigma[i], d, vector, fv);
if (gx > gmax)
{
predicted_class = i;
gmax = gx;
}
}
free(fv);
}
}
else
errno = EINVAL;
return predicted_class;
}
/*
Calculate quadratic classifier output for 'vector'. The classifier
is defined by 'model'. It is assumed that the length of 'vector' is
ad-1. The output is a vector of length c which contains values of
discriminant function for each class.
Return -1 in case of bad arguments or malloc() error and set errno.
*/
float *pqc_output(struct qmodel *model, int c, int ad, float *vector)
{
int i;
int d;
int predicted_class;
float gmax;
float *fv;
float *output;
output = (float *) 0;
predicted_class = -1;
gmax = -FLT_MAX;
d = ad-1;
if (model && (c > 0) && (ad > 0) && vector)
{
output = malloc(c*sizeof(float));
if (output)
{
fv = malloc(d*sizeof(float));
if (fv)
{
for (i = 0; i < c; i++)
output[i] = pqc_compute(model->wmx[i], model->sigma[i], d, vector, fv);
free(fv);
}
}
}
else
errno = EINVAL;
return output;
}
syntax highlighted by Code2HTML, v. 0.9.1