eval_svm.c

linux下的源码分类器SVM

/* Copyright 2001, 2002, 2003 Yann Guermeur and Andre Elisseeff               */

/* 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.                               */

/* You should have received a copy of the GNU General Public License          */
/* along with this program; if not, write to the Free Software                */
/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA  */

/*----------------------------------------------------------------------------*/
/*  Name           : eval_SVM.c                                               */
/*  Version        : 1.0                                                      */
/*  Creation       : 04/30/00                                                 */
/*  Last update    : 06/30/08                                                 */
/*  Subject        : Implementation of a M-SVM                                */
/*  Module         : M-SVM used in test                                       */
/*  Author         : Yann Guermeur [e-mail]                     */
/*----------------------------------------------------------------------------*/


#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include "algebre.h"
#include "biblio.h"

#define true 1
#define false 0
#define inf(a,b) ((a)<=(b)?(a):(b))
#define sup(a,b) ((a)<=(b)?(b):(a))

#define taille 81
#define tres_grand 1e6
#define tres_petit 1e-4
#define pas 10

FILE *fs, *fc;

char choice='y';

long nature_kernel=0, *y_app, *y_test,
Q=0, i, j, k, l, dim_input, nb_data_app=0, nb_data_test=0, y_i, y_j,
bad_b = false;

char fichier_app[taille], fichier_alpha[taille], fichier_SV[taille],
fichier_test[taille], fichier_hyp[taille], commande[taille],
conf[taille], fichier_fichcom[taille], fichier_resultat[taille],
*chaine_app =  "apprentissage",
*chaine_test = "test";

double **X_app, **X_test, **nb_sup_vect, R_emp=0.0,
**alpha, **delta_b_SVM, *b_SVM, C=0.0, **sup_vect, step=1e-1, partiel,
primal=0.0, dual=0.0, Q_form_dual=0.0, **gradient, negligeable, max_alpha=0.0,
**mat_conf, cost=0.0, *best_b_SVM, R_emp_min=0.0, check=0,
**H_alpha, **W, norm, **margin_vect, **nb_margin_vect, **radius;

/* Functions included in this program */

void caract_db();
void read_data();
void read_alpha();
void check_feasible_sol();
void compute_gradient();
void check_sup_vect();
void check_margin_vect();
double fonction(long categorie, double *vecteur);
void eval_training();
void compute_W();
void estime_b();
void compute_outputs(double **X, long *y, char *chaine, long nb_data);
void compute_obj_dual();
void compute_obj_primal();
void optimize_b();
void standardize_b();
void compute_boundaries();
void compute_significance();
void Pause(char *message);

main(int argc, char *argv[])

{

strcpy(fichier_fichcom, argv[1]);

system("clear");
caract_db();
read_data();
read_alpha();
check_feasible_sol();
compute_gradient();
check_sup_vect();
check_margin_vect();
if(nature_kernel == 1)
  compute_W();
estime_b();
eval_training();
/* printf("\nFirst estimate of the training performance"); */
compute_obj_dual();
compute_obj_primal();
if((max_alpha < C * (1.0 - tres_petit)) && (- Q_form_dual >= dual))
  {
  for(k=1; k<=Q; k++)
    check += mat_conf[k][k];
  if(check == nb_data_app)
    {
    bad_b = true;
    if(primal > -(1.5 * Q_form_dual))
      printf("\nBeware, very poor estimates of the bias (primal = %lf)\n",
      primal);
    primal = - Q_form_dual;
    }
  }

printf("\n*** Training performance\n");
printf("\nDifference between the primal and the dual: %lf -> %lf",
primal, dual);
printf("\n(ratio %5.2lf\%)\n", 100.0 * dual / primal);
display_stats(chaine_app, nb_data_app, nb_data_app, Q, mat_conf, cost);
display_full_stats(chaine_app, nb_data_app, Q, mat_conf, cost);
printf("Do you want vector b to be optimized [y/n]? ");
choice = getc(stdin);

if((choice == 'y') || (choice == 'Y'))
  {
  optimize_b();
  printf("\nDifference between the primal and the dual: %lf -> %lf",
  primal, dual);
  printf("\n(ratio %5.2lf\%)\n", 100.0 * dual / primal);
  if(bad_b == true)
    printf("\nWith optimal values for the biases, the primal should be: %lf\n",
    - Q_form_dual);
  display_stats(chaine_app, nb_data_app, nb_data_app, Q, mat_conf, cost);
  display_full_stats(chaine_app, nb_data_app, Q, mat_conf, cost);
  Pause("");
  }

if((nature_kernel == 1) && (dim_input == 2) && (Q == 3) && (max_alpha > 0.0))
  compute_boundaries();

printf("\n*** Test performance\n");
compute_outputs(X_test, y_test, chaine_test, nb_data_test);

display_full_stats(chaine_test, nb_data_test, Q, mat_conf, cost);
/* compute_significance(); */

}

void Pause(char *message)

{

char c;

static int n = 0;

if(message!=NULL)
  printf("%s\n",message);

printf("\a=> New Line to proceed...\n");

n = read(0, &c, 1);

while(c != '\n')
  n = read(0, &c, 1);

return;

}

void standardize_b()

{

double sum_b = 0.0;

for(k=1; k<=Q; k++)
  sum_b += b_SVM[k];

sum_b /= Q;

for(k=1; k<=Q; k++)
  b_SVM[k] -= sum_b;

}

void optimize_b()

{

long compteur, portion, class, Iteration;

double *grid_index, *b_SVM_base;

R_emp_min = R_emp;

for(class=1; class<=Q; class++)
  best_b_SVM[class] = b_SVM[class];

grid_index = (double *) calloc(Q+1, sizeof(double));
b_SVM_base = (double *) calloc(Q+1, sizeof(double));

grid_index[1] = 0.0;
b_SVM_base[1] = 0.0;

for(Iteration=1; Iteration <=10; Iteration++)
  {
  for(k=2; k<=Q; k++)
    b_SVM_base[k] = b_SVM[k] - (5.0 * step);

  for(compteur=0; compteur<(long)pow(10.0, (double)(Q-1)); compteur++)
    {
  
/*    printf("\nPosition in the grid: %6d\n", compteur); */

    portion = compteur;
    for(class=2; class<=Q; class++)
      {
      grid_index[class] = 
      (double) (portion / (long) pow(10.0, (double) (Q-class)));
      portion -= (long) (grid_index[class] * pow(10.0, (double) (Q-class)));
      }

    for(class=2; class<=Q; class++)
      b_SVM[class] = b_SVM_base[class] + (grid_index[class] * step);

    eval_training();
    compute_obj_primal();
    }

  for(class=1; class<=Q; class++)
    b_SVM[class] = best_b_SVM[class];

  eval_training();
  compute_obj_primal();

/*  
  scanf("%s", conf);
  sleep(2);
*/

  step /= 2.0;
  }

standardize_b();

printf("\nNew components of the b vector:\n");

printf("\n");
for(class=1; class<=Q; class++)
  printf("%11.6f\n", b_SVM[class]);

}

void caract_db()

{

if((fs=fopen(fichier_fichcom, "r"))==NULL)
  {
  printf("\nFile of parameters: %s cannot be open...\n", fichier_fichcom);
  exit(0);
  }

fscanf(fs, "%d", &Q);

mat_conf = matrix(Q,Q);

fscanf(fs, "%d", &nature_kernel);
fscanf(fs, "%lf", &C);

negligeable = C * tres_petit;

fscanf(fs, "%s", fichier_app);
printf("\nThe file of training is: %s\n", fichier_app);
fscanf(fs, "%s", fichier_test);
printf("\nThe file of test is: %s\n", fichier_test);
sleep(1);

fscanf(fs, "%s", fichier_alpha);
fscanf(fs, "%s", fichier_SV);
fscanf(fs, "%s", fichier_resultat);

fclose(fs);

}

void read_data()

{

long min_y = 1000;

if((fs=fopen(fichier_app, "r"))==NULL)
  {
  printf("\nFile of data: %s cannot be open...\n", fichier_app);
  exit(0);
  }

fscanf(fs, "%d", &nb_data_app);
fscanf(fs, "%d", &dim_input);

X_app = matrix(nb_data_app, dim_input);
y_app = (long *) calloc(nb_data_app+1, sizeof(long));

for(i=1; i<=nb_data_app; i++)
  {
  for(j=1; j<=dim_input; j++)
    fscanf(fs, "%lf", &X_app[i][j]);
  fscanf(fs, "%lf", &partiel);
  y_app[i] = (long) partiel;
  if(y_app[i] < min_y)
    min_y = y_app[i];
  }

fclose(fs);

if((fs=fopen(fichier_test, "r"))==NULL)
  {
  printf("\nFile of data: %s cannot be open...\n", fichier_test);
  exit(0);
  }

fscanf(fs, "%d", &nb_data_test);
fscanf(fs, "%d", &dim_input);

X_test = matrix(nb_data_test, dim_input);
y_test = (long *) calloc(nb_data_test+1, sizeof(long));

for(i=1; i<=nb_data_test; i++)
  {
  for(j=1; j<=dim_input; j++)
    fscanf(fs, "%lf", &X_test[i][j]);
  fscanf(fs, "%lf", &partiel);
  y_test[i] = (long) partiel;
  if(y_test[i] < min_y)
    min_y = y_test[i];
  }

fclose(fs);

if(min_y == 0)
  {
  printf("\nStandardizing the coding of the categories...\n");
  sleep(2);
  for(i=1; i<=nb_data_app; i++)
    y_app[i]++;
  for(i=1; i<=nb_data_test; i++)
    y_test[i]++;
  }

printf("\nThe files of training and test data have been read...\n");
sleep(1);

}

void read_alpha()

{

max_alpha = 0.0;

if((fs=fopen(fichier_alpha, "r"))==NULL)
  {
  printf("\nFile of dual variables: %s cannot be open...\n", 
  fichier_alpha);
  exit(0);
  }

alpha = matrix(nb_data_app, Q);

for(i=1; i<=nb_data_app; i++)
  for(k=1; k<=Q; k++)
    {
    fscanf(fs, "%lf", &alpha[i][k]);
    if((alpha[i][k] != 0.0) && (k == y_app[i]))
      {
      printf("\nEx. %d: pb. with the dummy variables...\n\n", i);
      exit(0);
      }
    if(alpha[i][k] > max_alpha)
      max_alpha = alpha[i][k];
    }

fclose(fs);

printf("\nThe maximum value among the dual variables is: %lf\n", max_alpha);
sleep(2);

}

void check_feasible_sol()

{

double *constraints;

for(i=1; i<=nb_data_app; i++)
  for(k=1; k<=Q; k++)
    if((alpha[i][k] < 0.0) || (alpha[i][k] > C))
       {
       printf("\nNo feasible solution: alpha[%d][%d] = %lf\n\n",
              i, k, alpha[i][k]);
       exit(0);
       }
     
constraints = calloc(Q+1, sizeof(double));

for(k=1; k<=Q; k++)
  constraints[k] = 0.0;

for(i=1; i<=nb_data_app; i++)
  for(k=1; k<=Q; k++)
     {
     if(y_app[i] == k)
       for(l=1; l<=Q; l++)
         constraints[k] -= alpha[i][l];
     else
       constraints[k] += alpha[i][k];  
     }

printf("\nSatisfaction of the equality constraints:\n\n");

for(k=1; k<=Q; k++)
  printf("%11.8f\n", constraints[k]);

norm = 0.0;

for(k=1; k<=Q; k++)
  norm += constraints[k] * constraints[k];

norm = sqrt(norm);

if(norm >= 1e-4)
  {
  printf("\nLarge deviation of the equality constraints...\n");
  exit(0);
  }

sprintf(commande, "cp %s Feasible/.", fichier_alpha);
system(commande);

}

void compute_gradient()

{

gradient = matrix(nb_data_app, Q);
H_alpha  = matrix(nb_data_app, Q);

printf("\nStart of the gradient computation...\n");

for(i=1; i<=nb_data_app; i++)
  {
  if((i%100) == 0)
    printf("\n%9d", i);
  y_i = y_app[i];

  for(k=1; k<=Q; k++)
    {
    H_alpha[i][k] = 0.0;
    if(k != y_i)
      {
      for(j=1; j<=nb_data_app; j++)