home *** CD-ROM | disk | FTP | other *** search
/ ftp.disi.unige.it / 2015-02-11.ftp.disi.unige.it.tar / ftp.disi.unige.it / pub / .person / BarlaA / sw / OLD / Simo / SVM_mono / solve.c < prev    next >
C/C++ Source or Header  |  2002-06-25  |  13KB  |  459 lines
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <time.h>
  5. #include <math.h>
  6. #include "patt_rec.h"
  7. #include "image2d.h"
  8. #include "pnmio.h"
  9. #include "nrutil.h"
  10.  
  11. #define x(i)        vect_training(train, (i))
  12. #define y(i)        class_training(train, (i))
  13. #define x_in(i)        vect_training(train, index_in[(i)])
  14. #define y_in(i)      class_training(train, index_in[(i)])
  15. #define x_out(i)      vect_training(train, index_out[(i)])
  16. #define y_out(i)     class_training(train, index_out[(i)])
  17. #define alpha_in(i)    alpha[index_in[(i)]]
  18. #define alpha_out(i)      alpha[index_out[(i)]]
  19. #define svc_in(i)         svc[index_in[(i)]]
  20. #define svc_out(i)        svc[index_out[(i)]]    
  21.  
  22. #define EPS             10e-6
  23.  
  24. static     int        *svc;        /* 0=non sv ; 1=margine ; 2=outlier    */
  25. static     int        *index_in;    /* Array degli indici dei vett. usati    */
  26. static    int          *index_out;    /* Array degli indici dei vett. out    */
  27. static     double        *alpha;        /* Array degli alpha        */
  28. static     training    *train;        /* Vettori di input usati da solve    */
  29.  
  30. extern int    knl;
  31. extern int    chunk_size;
  32. extern int    ell;
  33. extern int    dim;
  34. extern double    c_const;
  35. extern double    bee;
  36. extern int    degree;
  37. extern double    sigma, normalization;
  38. extern char     *image_ker;
  39. extern int      dim_r, dim_c;
  40. double          *x0;
  41. double          Radius;
  42. double          **mat_kernel;       /* MATRICE KERNEL PRECALCOLATO */
  43.  
  44. char            **mat2, **mat3;
  45. Image2D         *mat_image;
  46. int             maximum=-1, minimum=10;
  47. int             **mat;
  48. int             giro;
  49. FILE *fp, *f_ker;
  50.  
  51.  
  52. /* Se la soluzione e` ottimale ritorna 1, altrimenti ritorna 0. Per la verifica 
  53.    controlla (l - chunk_size) elementi del training set; se qualche 
  54.    punto viola le condizioni di ottimalita` (Khun-Tucker) aggiorna di 
  55.    conseguenza index_in ed index_out */
  56.    
  57. int optimal()
  58. {
  59.     register    i, j, margin_sv_number, error_sv_number, z, k, s;
  60.     int        nb = ell - chunk_size,    /* nb e` il numero di punti fuori base    */
  61.             nb_choose;        /* n.ro max. di p.ti out da testare    */
  62.  
  63.     double    gx_i,         /* parametro delle equazioni di KT    */
  64.             k_i,         /* altro parametro            */
  65.             aux, aux2, aux3, max, min;    
  66.  
  67.     double *point;
  68.  
  69.     nb_choose = ell - chunk_size;
  70.  
  71. //    fprintf(stderr,"dentro optimal\n");
  72.  
  73.     /* Computing Radius */
  74.     
  75.     aux = 0.;
  76.     margin_sv_number = 0;
  77.  
  78.     // calcolo il centro 
  79.       /*      for (i=0; i<dim; i++) 
  80.         x0[i] = 0.;     
  81.           
  82.       for (i = 0; i < ell; i++)
  83.         //if (svc_in(i) != 0) //if (alpha[i] > DELTA) 
  84.         {
  85.             point = x(i);
  86.             for (j = 0; j < dim; j++){
  87.             x0[j]+=alpha[i]*point[j];
  88.                   }
  89.         }      
  90.  
  91.     ker_c = kernel(knl,x0,x0);
  92.         */
  93.  
  94.     for (i = chunk_size - 1; i >= 0; i--)
  95.       if (svc_in(i) == 1) 
  96.         {
  97.           margin_sv_number++;
  98.               aux += mat_kernel[index_in[i]][index_in[i]];
  99.           for (j=ell-1;j>=0;j--) {
  100.         aux += -2.*alpha[j]*mat_kernel[index_in[i]][j]; 
  101.         for (k=ell-1;k>=0;k--)
  102.           aux += alpha[j]*alpha[k]*mat_kernel[j][k];
  103.           }
  104.         }
  105.  
  106.     //fprintf(stderr,"aux2 %lf\n",aux2);
  107.     
  108.     if (margin_sv_number == 0) 
  109.     {
  110.         error_sv_number = 0;
  111.         min = 2000000000.;
  112.         for (i = chunk_size - 1; i >= 0; i--) 
  113.           if (svc_in(i) == 2) 
  114.             {
  115.               aux= mat_kernel[index_in[i]][index_in[i]];
  116.               error_sv_number++;
  117.               for (j=ell-1;j>=0;j--) {
  118.             aux += -2.*alpha[j]*mat_kernel[index_in[i]][j];
  119.             for (k=ell-1;k>=0;k--)
  120.               aux += alpha[j]*alpha[k]*mat_kernel[j][k];
  121.               }
  122.  
  123.  
  124.               if (aux < min) 
  125.             min = aux;
  126.             }
  127.  
  128.         if (error_sv_number == 0) 
  129.         {
  130.             fprintf(stderr,"Strange: no support vector found!\n");
  131.             exit(1);
  132.         }
  133.         aux=min;
  134.         Radius = sqrt(min);        
  135.     }
  136.     else {
  137.     //  fprintf(stderr,"aux %lf\n",aux);
  138.       aux /= margin_sv_number;
  139.       Radius = sqrt(aux);
  140.     }
  141. //    fprintf(stderr,"dentro optimal...radius %lf\n",Radius);
  142.  
  143.     /* inserisco in base i punti che violano le condizioni di Kuhn-Tucker
  144.         scegliendo tra nb_choose non utilizzati (ne inserisco al massimo
  145.        chunk_size-2)  */
  146.     
  147.     srandom((unsigned) time(NULL));
  148.     for (z = nb_choose, s = chunk_size-2; z > 0 && s > 0; z--, nb--) 
  149.     {
  150.         /* pesco con probabilita` uniforme i punti fuori base */
  151.         i = random() % nb;
  152.         gx_i = aux - mat_kernel[index_out[i]][index_out[i]];  //aux raggio al quadrato
  153.  
  154.                          
  155.           for (j=ell-1;j>=0;j--) {
  156.         gx_i += 2.*alpha[j]*mat_kernel[index_out[i]][j];
  157.         for (k=ell-1;k>=0;k--)
  158.           gx_i -= alpha[j]*alpha[k]*mat_kernel[j][k];
  159.           }
  160.  
  161.                     
  162.         if (((alpha_out(i) < DELTA) && (gx_i < (-DELTA))) ||
  163.             ((alpha_out(i) > (c_const-DELTA)) && (gx_i > (+DELTA))) ||
  164.             ((alpha_out(i) > DELTA) && (alpha_out(i) < c_const-DELTA) && 
  165.              ((gx_i < -DELTA) || (gx_i > +DELTA)))) 
  166.         {
  167.             k = random() % s + 1;
  168.             j = index_in[k];
  169.             memmove(index_in+k, index_in+k+1, (s-k)*sizeof(int));
  170.             index_in[s--] = index_out[i];
  171.         } 
  172.         else
  173.             j = index_out[i];
  174.         memmove(index_out+i, index_out+i+1, (nb-i-1)*sizeof(int));
  175.         index_out[nb-1] = j;
  176.     }
  177.     //fprintf(stderr, "%d\n", chunk_size-s-2);
  178.     if (s == chunk_size-2)
  179.         return 1;
  180.     else
  181.         return 0;
  182. }
  183.  
  184. void  pre_compute(training *tr,int knl,double **mat_kernel)
  185. {
  186.         int i,j;
  187.  
  188.         train = tr;
  189.  
  190.         ell = length_training(train);
  191.  
  192.     
  193.         for(i=0;i<ell;i++){
  194.                 for(j=i;j<ell;j++){
  195.                         //printf("%d %d\n",i,j);
  196.                         mat_kernel[i][j] = kernel(knl, x(i), x(j)); 
  197.             fprintf(stderr,"%lf ",mat_kernel[i][j]);
  198.                        }
  199.           printf("\n");
  200.           }
  201.  
  202.         for(i=0;i<ell;i++)
  203.                 for(j=0;j<i;j++){
  204.                         mat_kernel[i][j] = mat_kernel[j][i];
  205.                         }
  206. }
  207.  
  208.  
  209. double *solve(training *tr, int knl, double c_const, int sp_dim, char *image_ker, char *f_kernel)
  210. {
  211.     register    i, j;
  212.     int         *sp_y,        /* Vettore delle classi         */
  213.             *sp_svc;    /* Classificazione dei vettori        */
  214.     double        **sp_D,     /* Componenti quadratiche del funz.    */
  215.             *sp_h,         /* Componenti lineari del funz.        */
  216.             *sp_alpha;
  217.     double         sp_e, aux, inc, min, max;
  218.     double         *solution;
  219.     
  220.      
  221.     train = tr;
  222. //    sp_dim *= 2;
  223.     ell = length_training(train);
  224.     dim = dim_training(train);
  225. //    chunk_size = (sp_dim+2 < ell) ? sp_dim+2 : ell;
  226.     chunk_size = (sp_dim+1 < ell) ? sp_dim+1 : ell;
  227.     fprintf(stderr,"chunksize %d ell %d\n",chunk_size,ell);
  228.  
  229.     /* Allocazione memoria necessaria ai dati    */
  230.     
  231.     index_in = (int *)myalloc(NULL, chunk_size*sizeof(int), "solve: index_in\n");
  232.     if (ell > chunk_size)
  233.         index_out = (int *)myalloc(NULL, (ell-chunk_size)*sizeof(int), "solve: index_out\n");
  234.     else
  235.         index_out = NULL;
  236.     alpha = (double *)myalloc(NULL, ell*sizeof(double), "solve: alpha\n");
  237.  
  238.     //XX x0 (dim)
  239.  
  240.     x0 = (double *)malloc(sizeof(double)*dim);
  241.   
  242.     svc = (int*)myalloc(NULL, ell*sizeof(int), "solve: svc\n");
  243.     for (i = ell-1; i >= 0; i--)
  244.         alpha[i] = 0.;
  245.     sp_alpha = (double *)myalloc(NULL, chunk_size*sizeof(double), "solve: sp_alpha\n");
  246.     sp_svc= (int *)myalloc(NULL, chunk_size*sizeof(int), "solve: sp_svc\n");
  247.     sp_D = (double **)myalloc(NULL, chunk_size*sizeof(double *), "solve: sp_D\n");
  248.     for (i = chunk_size - 1; i >= 0; i--)
  249.         sp_D[i]    = (double *)myalloc(NULL, chunk_size*sizeof(double), "solve: sp_D[i]\n");
  250.  
  251.     sp_h = (double *)myalloc(NULL, chunk_size*sizeof(double), "solve: sp_h\n");
  252.     sp_y = (int *) myalloc(NULL, chunk_size*sizeof(int), "solve: sp_cl\n");
  253.  
  254.     mat_kernel = dmatrix(0,ell-1,0,ell-1);
  255.  
  256.     // calcolo immagine solo se non faccio chunking
  257.     /*    
  258.       if (chunk_size==ell) {
  259.             mat = imatrix(0,ell-1,0,ell-1);
  260.             mat2 = cmatrix(0,ell-1,0,ell-1);
  261.         mat3 = cmatrix(0,ell-1,0,ell-1);
  262.       }
  263.     */
  264.            pre_compute(tr,knl,mat_kernel);
  265.  
  266.     f_ker = fopen(f_kernel,"w");
  267.     fprintf(f_ker,"%d \n",ell);
  268.     /*if (knl == NEGATIVE_KERNEL){
  269.       for(i=0;i<ell;i++) {
  270.         for(j=0;j<ell;j++)
  271.           fprintf(f_ker,"%lf ",-mat_kernel[i][j]);
  272.         fprintf(f_ker,"\n");
  273.       }
  274.     }
  275.     else{*/
  276.       for(i=0;i<ell;i++) {
  277.         for(j=0;j<ell;j++)
  278.           fprintf(f_ker,"%lf ",mat_kernel[i][j]);
  279.         fprintf(f_ker,"\n");
  280.       }
  281. //    }
  282.     fclose(f_ker);
  283.  
  284.  
  285.     /* Creazione del sottoproblema: 
  286.          scelta dei primi indici da utilizzare e da non utilizzare    */
  287.  
  288.         // modifica: noi abbiamo solo una classe....
  289.  
  290.     for (i = 0, j = 0; i < chunk_size; i++, j++)
  291.               index_in[j] = i;          /* I primi chunk_size di classe 1     */
  292.     
  293.     for (i = chunk_size, j = 0; i < ell; i++, j++)
  294.         index_out[j] = i;
  295.  
  296.     /* Inizio del ciclo di risoluzione del problema    */
  297.  
  298.         fprintf(stderr,"prima del do\n");
  299.     giro=0;
  300.     do 
  301.     {
  302.       giro++;
  303.       fprintf(stderr,"giro %d\n",giro);
  304.       // debug
  305.  
  306.     //    fprintf(stderr,"index_in\n");
  307.      // for (i=0;i<chunk_size;i++)
  308.        // fprintf(stderr,"%d ",index_in[i]);
  309.       //fprintf(stderr,"\n");
  310.       
  311.       //fprintf(stderr,"index_out\n");
  312.       //for (i=0;i<ell-chunk_size;i++)
  313.         //fprintf(stderr,"%d ",index_out[i]);
  314.       //fprintf(stderr,"\n");
  315.       
  316.       /* Preparazione delle variabili da passare al solver    */
  317.       
  318.       /* Matrice della parte quadratica della funzione obbiettivo    */
  319.       
  320.       for (i = chunk_size-1; i >= 0; i--){
  321.         //      fprintf(stderr,"%d\n",i);
  322.         for (j = chunk_size-1; j >= 0; j--) {
  323.           sp_D[i][j] =  y_in(i)*y_in(j)*mat_kernel[index_in[i]][index_in[j]];
  324.             
  325.         }
  326.       }
  327.       /* if (chunk_size == ell) {
  328.          minimum = -maximum;
  329.          mat_image = ImageAlloc(ell,ell);
  330.          for (i = chunk_size-1; i >= 0; i--)
  331.          for (j = chunk_size-1; j >= 0; j--) {
  332.          mat2[i][j] = ( mat[i][j] - minimum )*255/(maximum - minimum);
  333.          mat_image->data[i][j]=mat2[i][j];
  334.          }         
  335.          
  336.          write_file_pnm(image_ker, mat_image );
  337.          }*/
  338.       
  339.       //fprintf(stderr,"prima del calcolo lineare\n");
  340.       
  341.       /* Vettore della parte lineare della funzione obbiettivo    */
  342.       for (i = chunk_size-1; i >= 0; i--) 
  343.         {
  344.           aux = 0.;
  345.           for (j = ell-chunk_size-1; j >= 0; j--)
  346.         aux += 2*y_out(j)*alpha_out(j)*mat_kernel[index_in[i]][index_out[j]]; //attenti al 2!
  347.           //sp_h[i] = y_in(i) * aux - 1.;
  348.           sp_h[i] = y_in(i) * aux - mat_kernel[index_in[i]][index_in[i]];
  349.         }
  350.       
  351.       //fprintf(stderr,"prima dei vincoli\n");
  352.       
  353.       /* Costante del primo gruppo di vincoli    */
  354.       //sp_e = 0.;
  355.       sp_e = -1.;
  356.       //fprintf(stderr,"alpha fuori\n");
  357.       for (i = ell-chunk_size-1; i >= 0; i--) {
  358.         sp_e += y_out(i) * alpha_out(i); //mat_kernel[index_out[i]][index_out[i]] * alpha_out(i);
  359.         //fprintf(stderr,"%lf ",alpha_out(i));
  360.       }
  361.       
  362. //      fprintf(stderr,"\n funzionavano meglio le parolacce. sp_e %lf\n",sp_e);
  363.       
  364.     /* Vettore delle classi dei punti in base    */
  365.         for (i = chunk_size-1; i >= 0; i--)
  366.             //sp_y[i] = mat_kernel[index_in[i]][index_in[i]];
  367.       sp_y[i] = y_in(i);
  368.       
  369.       /* Chiamate al solver */
  370.       
  371.       fprintf(stderr,"prima del solver\n");
  372.       solver(chunk_size, sp_D, sp_h, c_const, sp_e, sp_y, sp_alpha, sp_svc);
  373.       
  374.       fprintf(stderr,"dopo il solver\n");
  375.       
  376.       /* Aggiornamento del vettore globale delle soluzioni */
  377.       
  378.     //  fprintf(stderr,"alpha in\n");    
  379.       for (i = chunk_size-1; i >= 0; i--) 
  380.         {
  381.           alpha_in(i) = sp_alpha[i];
  382. //          fprintf(stderr,"%lf ",alpha_in(i));
  383.           svc_in(i) = sp_svc[i];
  384.         }
  385.       fprintf(stderr,"\n");
  386.       
  387.     } while (!optimal());
  388.     
  389.     solution = alpha;
  390.     
  391.     
  392.         free(index_in);
  393.     free(index_out);
  394.     for (i = chunk_size - 1; i >= 0; i--)
  395.       free(sp_D[i]);
  396.      free(sp_D);
  397.      free(sp_alpha);
  398.     free(sp_svc);
  399.     
  400.     return solution;
  401.       }
  402.  
  403. void write_solution(FILE *fp, double *sol, int knl, double bee, double c_const)
  404. {
  405.   register    i, j;
  406.   int        n,m,temp,sv_number = 0;
  407.   double    alpha_sum, *point;
  408.   double          *R,*quad,R_alpha,R_alpha2,R_sum;
  409.   
  410.   alpha_sum=0.;
  411.   
  412.   for (i = 0; i < ell; i++){
  413.     if (alpha[i] > DELTA) sv_number++;
  414.     alpha_sum+=alpha[i];
  415.   }
  416.   
  417.   fprintf(stderr, "sv_number: %d\n", sv_number);
  418.   fprintf(fp, "%d ", sv_number);
  419.   fprintf(stderr,"dimensione punti: %d\n", dim);
  420.   fprintf(fp, "%d ", dim);
  421.   fprintf(fp, "%d %d\n",dim_r,dim_c);
  422.   
  423.   
  424.   //Scrivo gli sv, le labels e gli alpha e x0
  425.     
  426.     
  427.     for (i = 0; i < ell; i++)
  428.       // fprintf(stderr,"i: %d alpha: %lf\n",i,alpha[i]);
  429.   if (svc[i]!=0) 
  430.     {
  431.       point = x(i);
  432.       //modifica al file dei sv: aggiungiamo gli indici
  433.       fprintf(fp,"%d ",i);
  434.       for (j = 0; j < dim; j++){
  435.     fprintf(fp, "%lf ", point[j]);
  436.       }
  437.       fprintf(fp, "%d ", y(i));
  438.       fprintf(fp, "%lf\n", alpha[i]);
  439.       printf("i: %d alpha: %lf\n",i,alpha[i]);
  440.     }             
  441.     
  442.   fprintf(fp, "%lf\n",Radius);
  443.   
  444.   
  445. //Ora la costante C, il tipo di kernel e la normalizzazione
  446.   fprintf(fp, "%lf\n", c_const);    
  447.   fprintf(fp, "%d\n", knl);
  448.   if (knl == RBF_KERNEL)
  449.     fprintf(fp, "%lf\n", sigma);
  450.   else if (knl == POLY_KERNEL)
  451.     fprintf(fp, "%lf %d\n", normalization, degree); 
  452.   else if (knl == HAUS_KERNEL)
  453.     fprintf(fp, "%lf\n", normalization);
  454.   else 
  455.     fprintf(fp, "%lf\n", normalization);
  456.     
  457. }
  458.  
  459.