svm_matrix.cpp

Go to the documentation of this file.

#include <vector>
#include <string>
#include <fstream>
#include <iostream>

/* uncomment, if you want to use svm-learn out of C++ */
extern "C" { 
# include "svm_common.h"
# include "plugin.h"
}
#include <limits.h>

#define fixed_point_shift USHRT_MAX
#define fixed_point_type unsigned short

using namespace std;

void   read_input_parameters(int, char **, char *, char *, char *,
                             long *, long *, KERNEL_PARM *,
                             long *, long *, long *);
                             
void   wait_any_key();
void   print_help();

void setValue(void * matrix, long i, long j, float value, long no_cols, long fixed_point)
{
  if(!no_cols) // symmetric
    if(fixed_point)
      ((fixed_point_type*)matrix)[j+i*(i+1)/2] = (fixed_point_type)(value * fixed_point_shift);  
    else
      ((float*)matrix)[j+i*(i+1)/2] = value;
  else // asymmetric
    if(fixed_point)
      ((fixed_point_type*)matrix)[j+i*no_cols] = (fixed_point_type)(value * fixed_point_shift);
    else
      ((float*)matrix)[j+i*no_cols] = value;
}

float getValue(void * matrix, long i, long j, long no_cols, long fixed_point)
{
  if(!no_cols) // symmetric
    if(fixed_point)
      return ((float*)matrix)[j+i*(i+1)/2]/fixed_point_shift;
    else
      return ((float*)matrix)[j+i*(i+1)/2];
  else // asymmetric
    if(fixed_point)
      return ((float*)matrix)[j+i*no_cols]/fixed_point_shift;
    else
      return ((float*)matrix)[j+i*no_cols];
}

void printValue(FILE *fp, 
                void *value,
                long offset,
                long fixed_point, 
                long binary)
{  
  if(binary)
    if(fixed_point)
      fwrite(((fixed_point_type*)value)+offset, sizeof(fixed_point_type), 1, fp);
    else
      fwrite(((float*)value)+offset, sizeof(float), 1, fp);
  else
    if(fixed_point)
      fprintf(fp, "%hd ", *(((fixed_point_type*)value)+offset));
    else
      fprintf(fp, "%f ", *((float*)value+offset));
}

/* void * symmetric_matrix(char * docfile,  */
/*                      long verbosity,  */
/*                      long normalize,  */
/*                      KERNEL_PARM * kernel_parm_p, */
/*                      long fixed_point, */
/*                      long* no_rows) */
/* { */
/*   DOC **docs;  /\* training examples *\/ */
/*   long totwords,totdoc,i,j; */
/*   double *target; */
/*   void *matrix; */

/*   read_documents(docfile,&docs,&target,&totwords,&totdoc); */

/*   if(fixed_point) */
/*     matrix = malloc(sizeof(fixed_point_type)*totdoc*(totdoc+1)/2); */
/*   else */
/*     matrix = malloc(sizeof(float)*totdoc*(totdoc+1)/2); */
  
  
/*   if(normalize){ */
/*     float * diag = malloc(sizeof(float)*totdoc); */
/*     for(i = 0; i < totdoc; i++) */
/*       diag[i]=kernel(kernel_parm_p, docs[i], docs[i]); */
/*     for(i = 0; i < totdoc; i++) */
/*       for(j = 0; j < i; j++){ */
/*      float val = kernel(kernel_parm_p, docs[i], docs[j]); */
/*      float den = sqrt(diag[i]*diag[j]); */
/*      if(den != 0) */
/*        val /= den; */
/*      setValue(matrix,i,j, val, 0, fixed_point); */
/*       } */
/*     for(i = 0; i < totdoc; i++) // set diagonal elements to one */
/*       setValue(matrix,i,i,1,0,fixed_point);     */
/*     free(diag); */
/*   } */
/*   else */
/*     for(i = 0; i < totdoc; i++) */
/*       for(j = 0; j <= i; j++) */
/*      setValue(matrix,i,j, kernel(kernel_parm_p, docs[i], docs[j]), 0, fixed_point); */

/*   for(i=0;i<totdoc;i++)  */
/*     free_example(docs[i],1); */
/*   free(docs); */
/*   free(target); */

/*   *no_rows = totdoc; */
/*   return matrix; */
/* } */

/* void * asymmetric_matrix(char * rowfile, */
/*                       char * colfile, */
/*                       long verbosity,  */
/*                       long normalize,  */
/*                       KERNEL_PARM * kernel_parm_p, */
/*                       long fixed_point, */
/*                       long *no_rows, */
/*                       long *no_cols) */
/* { */
/*   DOC **row_docs,**col_docs;  /\* examples *\/ */
/*   long row_totwords,col_totwords,row_totdoc,col_totdoc,i,j; */
/*   double *row_target,*col_target; */
/*   void *matrix; */

/*   read_documents(rowfile,&row_docs,&row_target,&row_totwords,&row_totdoc); */
/*   read_documents(colfile,&col_docs,&col_target,&col_totwords,&col_totdoc); */


/*   if(fixed_point) */
/*     matrix = malloc(sizeof(fixed_point_type)*row_totdoc*col_totdoc); */
/*   else */
/*     matrix = malloc(sizeof(float)*row_totdoc*col_totdoc); */

  
/*   if(normalize){ */
/*     float * row_diag = malloc(sizeof(float)*row_totdoc); */
/*     float * col_diag = malloc(sizeof(float)*col_totdoc); */
/*     for(i = 0; i < row_totdoc; i++) */
/*       row_diag[i]=kernel(kernel_parm_p, row_docs[i], row_docs[i]); */
/*     for(i = 0; i < col_totdoc; i++) */
/*       col_diag[i]=kernel(kernel_parm_p, col_docs[i], col_docs[i]); */
/*     for(i = 0; i < row_totdoc; i++) */
/*       for(j = 0; j < col_totdoc; j++){ */
/*      float val = kernel(kernel_parm_p, row_docs[i], col_docs[j]); */
/*      float den = sqrt(row_diag[i]*col_diag[j]); */
/*      if(den != 0) */
/*        val /= den; */
/*      setValue(matrix, i, j, val, col_totdoc, fixed_point); */
/*       } */
/*     free(row_diag); */
/*     free(col_diag);     */
/*   } */
  
/*   for(i=0;i<row_totdoc;i++)  */
/*     free_example(row_docs[i],1); */
/*   free(row_docs); */
/*   free(row_target); */
/*   for(i=0;i<col_totdoc;i++)  */
/*     free_example(col_docs[i],1); */
/*   free(col_docs); */
/*   free(col_target); */

/*   *no_rows = row_totdoc; */
/*   *no_cols = col_totdoc; */
/*   return matrix; */
/* } */

void free_documents(DOC ** docs, double * target, long totdoc)
{
  for(long i=0;i<totdoc;i++) 
    free_example(docs[i],1);
  free(docs);
  free(target);
}

long getNoRows(const vector<string>& docfile)
{
  long norows = 0;
  ifstream in;
  string buf;
  
  for(uint i = 0; i < docfile.size(); i++){
    in.open(docfile[i].c_str());
    while(getline(in,buf))
      norows++;
    in.close();
  }
  return norows;
}

float * symmetric_matrix(void * matrix, DOC ** docs, long totdoc, long offset,
                         KERNEL_PARM * kernel_parm_p, long normalize, long fixed_point)
{
  float * diag = NULL;

  if(normalize){
    diag = (float*)malloc(sizeof(float)*totdoc);
    for(long i = 0; i < totdoc; i++)
      diag[i]=kernel(kernel_parm_p, docs[i], docs[i]);
    for(long i = 0; i < totdoc; i++)
      for(long j = 0; j < i; j++){
        float val = kernel(kernel_parm_p, docs[i], docs[j]);
            float den = sqrt(diag[i]*diag[j]);
            if(den != 0)
              val /= den;
            setValue(matrix,i+offset,j+offset, val, 0, fixed_point);
      }
    for(long i = 0; i < totdoc; i++) // set diagonal elements to one
      setValue(matrix,i+offset,i+offset,1,0,fixed_point);    
  }
  else
    for(long i = 0; i < totdoc; i++)
      for(long j = 0; j <= i; j++)
        setValue(matrix,i+offset,j+offset, kernel(kernel_parm_p, docs[i], docs[j]), 0, fixed_point);
  
  return diag;
}

void symmetric_matrix(void * matrix, float * diag_c, DOC ** docs_r, long totdoc_r, long offset_r,
                      DOC ** docs_c, long totdoc_c, long offset_c, 
                      KERNEL_PARM * kernel_parm_p, long normalize, long fixed_point)
{
  if(normalize){
    float * diag_r = (float*)malloc(sizeof(float)*totdoc_r);
    for(long i = 0; i < totdoc_r; i++)
      diag_r[i]=kernel(kernel_parm_p, docs_r[i], docs_r[i]);
    for(long i = 0; i < totdoc_r; i++)
      for(long j = 0; j < totdoc_c; j++){
        float val = kernel(kernel_parm_p, docs_r[i], docs_c[j]);
        float den = sqrt(diag_r[i]*diag_c[j]);
        if(den != 0)
          val /= den;
        setValue(matrix,i+offset_r,j+offset_c, val, 0, fixed_point);
      }
    free(diag_r);
  }
  else
    for(long i = 0; i < totdoc_r; i++)
      for(long j = 0; j < totdoc_c; j++)
        setValue(matrix,i+offset_r,j+offset_c, kernel(kernel_parm_p, docs_r[i], docs_c[j]), 0, fixed_point);  
}

void * symmetric_matrix(const vector<string>& docfile, 
                        long verbosity, 
                        long normalize, 
                        KERNEL_PARM * kernel_parm_p,
                        long fixed_point,
                        long* no_rows)
{
  uint blocks = docfile.size();
  DOC **docs[blocks];  /* training examples */
  long totwords[blocks],totdoc[blocks],fulltotdoc,offset_r = 0, offset_c = 0;
  double *target[blocks];
  void *matrix;
  
  fulltotdoc = getNoRows(docfile);

  if(fixed_point)
    matrix = malloc(sizeof(fixed_point_type)*fulltotdoc*(fulltotdoc+1)/2);
  else
    matrix = malloc(sizeof(float)*fulltotdoc*(fulltotdoc+1)/2);
 
  for(uint c = 0; c < docfile.size(); c++){
    read_documents(const_cast<char *>(docfile[c].c_str()),&docs[c],&target[c],&totwords[c],&totdoc[c]);  
    float * diag_c = symmetric_matrix(matrix,docs[c],totdoc[c],offset_c, kernel_parm_p, normalize, fixed_point);      
    offset_r = offset_c + totdoc[c];
    for(uint r = c+1; r < docfile.size(); r++){
      read_documents(const_cast<char *>(docfile[r].c_str()),&docs[r],&target[r],&totwords[r],&totdoc[r]);  
      symmetric_matrix(matrix, diag_c, docs[r], totdoc[r], offset_r, docs[c], totdoc[c], offset_c, 
                       kernel_parm_p, normalize, fixed_point);
      free_documents(docs[r],target[r],totdoc[r]);
      offset_r += totdoc[r];
    }
    free_documents(docs[c],target[c],totdoc[c]);
    if(diag_c) free(diag_c);
    offset_c += totdoc[c];
  }

  *no_rows = fulltotdoc;
  return matrix;
}

void asymmetric_matrix(void * matrix, DOC ** docs_r, long totdoc_r, long offset_r,
                       DOC ** docs_c, long totdoc_c, long offset_c, long nocols, 
                       KERNEL_PARM * kernel_parm_p, long normalize, long fixed_point)
{
  if(normalize){
    float * diag_r = (float*)malloc(sizeof(float)*totdoc_r);
    float * diag_c = (float*)malloc(sizeof(float)*totdoc_c);
    for(long i = 0; i < totdoc_r; i++)
      diag_r[i]=kernel(kernel_parm_p, docs_r[i], docs_r[i]);
    for(long i = 0; i < totdoc_c; i++)
      diag_c[i]=kernel(kernel_parm_p, docs_c[i], docs_c[i]);
    for(long i = 0; i < totdoc_r; i++)
      for(long j = 0; j < totdoc_c; j++){
        float val = kernel(kernel_parm_p, docs_r[i], docs_c[j]);
        float den = sqrt(diag_r[i]*diag_c[j]);
        if(den != 0)
          val /= den;
        setValue(matrix, i+offset_r, j+offset_c, val, nocols, fixed_point);
      }
    free(diag_r);
    free(diag_c);
  }
  else
    for(long i = 0; i < totdoc_r; i++)
      for(long j = 0; j < totdoc_c; j++)
        setValue(matrix,i+offset_r,j+offset_c, kernel(kernel_parm_p, docs_r[i], docs_c[j]), nocols, fixed_point);
}

void * asymmetric_matrix(const vector<string>& docfile_r,
                         const vector<string>& docfile_c,
                         long verbosity, 
                         long normalize, 
                         KERNEL_PARM * kernel_parm_p,
                         long fixed_point,
                         long *no_rows,
                         long *no_cols)
{
  uint blocks_r = docfile_r.size();
  DOC **docs_r[blocks_r];  /* training examples */
  long totwords_r[blocks_r],totdoc_r[blocks_r],fulltotdoc_r,offset_r = 0;
  double *target_r[blocks_r];
  uint blocks_c = docfile_c.size();
  DOC **docs_c[blocks_c];  /* training examples */
  long totwords_c[blocks_c],totdoc_c[blocks_c],fulltotdoc_c;
  double *target_c[blocks_c];
  void *matrix;
  
  fulltotdoc_r = getNoRows(docfile_r);
  fulltotdoc_c = getNoRows(docfile_c);

  if(fixed_point)
    matrix = malloc(sizeof(fixed_point_type)*fulltotdoc_r*fulltotdoc_c);
  else
    matrix = malloc(sizeof(float)*fulltotdoc_r*fulltotdoc_c);

  for(uint r = 0; r < blocks_r; r++){
    read_documents(const_cast<char *>(docfile_r[r].c_str()),&docs_r[r],&target_r[r],&totwords_r[r],&totdoc_r[r]);  
    long offset_c = 0;
    for(uint c = 0; c < blocks_c; c++){
      read_documents(const_cast<char *>(docfile_c[c].c_str()),&docs_c[c],&target_c[c],&totwords_c[c],&totdoc_c[c]);  
      asymmetric_matrix(matrix, docs_r[r], totdoc_r[r], offset_r, docs_c[c], totdoc_c[c], offset_c, 
                        fulltotdoc_c, kernel_parm_p, normalize, fixed_point);
      free_documents(docs_c[c],target_c[c],totdoc_c[c]);
      offset_c += totdoc_c[c];
    }
    free_documents(docs_r[r],target_r[r],totdoc_r[r]);
    offset_r += totdoc_r[r];
  }

  *no_rows = fulltotdoc_r;
  *no_cols = fulltotdoc_c;
  return matrix;
}


void print_matrix(const char * matrixfile, void * matrix, long fixed_point, 
                  long no_rows, long no_cols, long binary, long lower_traingular)
{
  FILE* fp;
  int i,j;

  fp = fopen(matrixfile,"w");
  if(fp == NULL){
    perror("ERROR in opening matrix file:");
    exit(1);
  }

  if(!no_cols) // symmetric
    if(lower_traingular)
      for(i = 0; i < no_rows; i++){
        for(j = 0; j <= i; j++)
          printValue(fp, matrix, (j+i*(i+1)/2), fixed_point, binary);
        if(!binary)
          fprintf(fp,"\n");
      }
    else
      for(i = 0; i < no_rows; i++){
        for(j = 0; j < no_rows; j++)
          printValue(fp, matrix, (j<=i) ? (j+i*(i+1)/2) : (i+j*(j+1)/2), fixed_point, binary);
        if(!binary)
          fprintf(fp,"\n");
      }
  else // asymmetric
    for(i = 0; i < no_rows; i++){
      for(j = 0; j < no_cols; j++)
        printValue(fp, matrix, (j+i*no_cols), fixed_point, binary);
      if(!binary)
        fprintf(fp,"\n");
    }    
  
  fclose(fp);
}

void getStringList(string buf, vector<string>& list)
{
  int beg = 0, end = 0;
  while((end = buf.find(":",beg)) != string::npos){
    list.push_back(buf.substr(beg,end-beg));
    beg = end+1;
  }
  list.push_back(buf.substr(beg));
}

/*---------------------------------------------------------------------------*/

void read_input_parameters(int argc,char *argv[],string& matrixfile, vector<string>& docfile_r,
                           vector<string>& docfile_c, long *verbosity,long *normalize,
                           KERNEL_PARM *kernel_parm, long *fixed_point, long *binary, 
                           long *lower_triangular)
{
  long i;
  
  /* set default */  
  strcpy (kernel_parm->plugin_file, "");
  (*verbosity)=1;
  (*normalize)=0;
  kernel_parm->kernel_type=0;
  kernel_parm->poly_degree=3;
  kernel_parm->rbf_gamma=1.0;
  kernel_parm->coef_lin=1;
  kernel_parm->coef_const=1;
  strcpy(kernel_parm->custom,"empty");
  (*fixed_point) = 0;
  (*binary) = 0;
  (*lower_triangular) = 0;
  
  for(i=1;(i<argc) && ((argv[i])[0] == '-');i++) {
    switch ((argv[i])[1]) 
      { 
      case '?': print_help(); exit(0);
      case 'v': i++; (*verbosity)=atol(argv[i]); break;
      case 'R': i++; getStringList(argv[i], docfile_r); break;
      case 'C': i++; getStringList(argv[i], docfile_c); break;
      case 't': i++; kernel_parm->kernel_type=atol(argv[i]); break;
      case 'd': i++; kernel_parm->poly_degree=atol(argv[i]); break;
      case 'g': i++; kernel_parm->rbf_gamma=atof(argv[i]); break;
      case 's': i++; kernel_parm->coef_lin=atof(argv[i]); break;
      case 'r': i++; kernel_parm->coef_const=atof(argv[i]); break;
      case 'u': i++; strcpy(kernel_parm->custom,argv[i]); break;
      case 'D': i++; strcpy(kernel_parm->plugin_file,argv[i]); 
        if (kernel_parm->kernel_type==0) kernel_parm->kernel_type=4;
        break;
      case 'N': (*normalize)=1; break;
      case 'B': (*binary)=1; break;
      case 'F': (*fixed_point)=1; break;        
      case 'L': (*lower_triangular)=1; break;   
      default: printf("\nUnrecognized option %s!\n\n",argv[i]);
               print_help();
               exit(0);
      }
  }
  if(i>=argc) {
    printf("\nNot enough input parameters!\n\n");
    wait_any_key();
    print_help();
    exit(0);
  }
  if(*fixed_point && !(*normalize)){
    printf("\nFixed point can only be used with normalized kernels!\n\n");
    wait_any_key();
    print_help();
    exit(0);
  }
  matrixfile = argv[i++];
}

void wait_any_key()
{
  printf("\n(more)\n");
  (void)getc(stdin);
}

void print_help()
{
  printf("\nSVM-light %s: Support Vector Machine, Gram matrix computation module     %s\n",VERSION,VERSION_DATE);
  copyright_notice();
  printf("   usage: svm_Dmatrix [options] matrix_file\n\n");
  printf("Arguments:\n");
  printf("         matrix_file-> file where kernel matrix will be saved\n");
  printf("General options:\n");
  printf("         -?                -> this help\n");
  printf("         -v [0..3]         -> verbosity level (default 1)\n");
  printf("         -R file1:..:fileN -> example files for rows (and cols if no -C -> i.e. symmetric)\n");
  printf("         -C file1:..:fileN -> example files for cols (asymmetric matrix)\n");
  printf("Kernel options:\n");
  printf("         -t int      -> type of kernel function:\n");
  printf("                        0: linear (default if -D not given)\n");
  printf("                        1: polynomial (s a*b+c)^d\n");
  printf("                        2: radial basis function exp(-gamma ||a-b||^2)\n");
  printf("                        3: sigmoid tanh(s a*b + c)\n");
  printf("                        4: user defined kernel as specified in plugin (default when -D is given)\n");
  printf("         -d int      -> parameter d in polynomial kernel\n");
  printf("         -g float    -> parameter gamma in rbf kernel\n");
  printf("         -s float    -> parameter s in sigmoid/poly kernel\n");
  printf("         -r float    -> parameter c in sigmoid/poly kernel\n");
  printf("         -u string   -> parameter of user defined kernel\n");
  printf("         -N          -> normalize kernel (default=false)\n");
  printf("         -D string   -> path to a dynamic library containing code for\n");
  printf("                        custom data types and associated kernel\n");
  printf("         -B          -> output matrix in binary form (instead of ascii)\n");
  printf("         -F          -> use fixed point values (needs normalized kernel)\n");
  printf("         -L          -> output lower triangular matrix only\n");
  wait_any_key();
  printf("\nMore details in:\n");
  printf("[1] T. Joachims, Making Large-Scale SVM Learning Practical. Advances in\n");
  printf("    Kernel Methods - Support Vector Learning, B. Schölkopf and C. Burges and\n");
  printf("    A. Smola (ed.), MIT Press, 1999.\n");
  printf("[2] T. Joachims, Estimating the Generalization performance of an SVM\n");
  printf("    Efficiently. International Conference on Machine Learning (ICML), 2000.\n");
  printf("[3] T. Joachims, Transductive Inference for Text Classification using Support\n");
  printf("    Vector Machines. International Conference on Machine Learning (ICML),\n");
  printf("    1999.\n");
  printf("[4] K. Morik, P. Brockhausen, and T. Joachims, Combining statistical learning\n");
  printf("    with a knowledge-based approach - A case study in intensive care  \n");
  printf("    monitoring. International Conference on Machine Learning (ICML), 1999.\n");
  printf("[5] T. Joachims, Learning to Classify Text Using Support Vector\n");
  printf("    Machines: Methods, Theory, and Algorithms. Dissertation, Kluwer,\n");
  printf("    2002.\n\n");
}

int main (int argc, char* argv[])
{  
  DOC **docs;  /* training examples */
  long totwords,totdoc,i,j;
  double *target;
  KERNEL_PARM kernel_parm;
  long normalize;
  long fixed_point;
  long binary;
  long lower_triangular;
  void *matrix;
  long no_rows = 0, no_cols = 0;
  string matrixfile;
  vector<string> docfile_r;
  vector<string> docfile_c;

  read_input_parameters(argc,argv,matrixfile,docfile_r,docfile_c,&verbosity,&normalize,&kernel_parm, 
                        &fixed_point, &binary, &lower_triangular);
  plugin_init(kernel_parm.plugin_file);
  plugin_kernel_setparm(kernel_parm.custom);
  

  if(docfile_c.size() == 0)
    matrix = symmetric_matrix(docfile_r,verbosity,normalize,&kernel_parm,fixed_point, &no_rows); 
  else
    matrix = asymmetric_matrix(docfile_r,docfile_c,verbosity,normalize,&kernel_parm,fixed_point, &no_rows, &no_cols); 

  print_matrix(matrixfile.c_str(), matrix, fixed_point, no_rows, no_cols, binary, lower_triangular);
  free(matrix);
  return(0);
}

---