algebre.c

linux下的源码分类器SVM

/*----------------------------------------------------------------------*/
/*  Version        : 1.0                                                */
/*  Creation       : 22/12/95                                           */
/*  Derniere modif.: 02/09/01                                           */
/*  Sujet          : Unites contenant des procedures d'algebre          */
/*  Intitule       : Librairie mathematique supplementaire II           */
/*  Auteur         : Yann Guermeur                                      */
/*----------------------------------------------------------------------*/

#include <stdio.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>

#define true 1
#define false 0
#define sign(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
#define NR_END 1
#define FREE_ARG char*
#define small 1e-8

char conf[81];

int i , j , k , l , dim , dim1 , dim2 , dim3;

double **mat1 , **mat2 , **mat3 , partiel , drand48();


void free_matrix(double **m)

/* Libere l'espace memoire occupe par une matrice */

{

free((double*)m[1]);
free((double*)m);

}

void dir_prod(double **mat1 , double **mat2 , int dim)

/* Calcul de vvT ou v est un vecteur colonne */

{

for(i=1;i<=dim;i++)
  for(j=1;j<=dim;j++)
    mat2[i][j] = mat1[i][1] * mat1[j][1];

}

void init_pos(double **mat , int dim1 , int dim2)

/* Initialise positivement les composantes d'une matrice */

{

srand48(getpid());

for(i=1;i<=dim1;i++)
  for(j=1;j<=dim2;j++)
    mat[i][j] = 0.1 * drand48() + 0.05;

}

double pythag(double a , double b)

/* Cf. Numerical Recipes in C p. 70 */

{

double result;

result = sqrt(a*a + b*b);

return result;

}

void tqli(double d[] , double e[], int dim , double **mat)

/* Valeurs propres des matrices symetriques tridiagonales */

{

double pythag(double a, double b);
int m,l,iter,i,k;
double s,r,p,g,f,dd,c,b;

for(i=2;i<=dim;i++)
  e[i-1] = e[i];

e[dim] = 0.0;

for(l=1;l<=dim;l++)
  {
  iter=0;
  do
    {
    for(m=l;m<=dim-1;m++)
      {
      dd = fabs(d[m]) + fabs(d[m+1]);
      if((double)(fabs(e[m])+dd) == dd)
        break;
      }
    if(m != l)
      {
      if(iter++ == 1000)
        {
        printf("\nToo many iterations in tqli...\n\n");
        exit(0);
        }
      g = (d[l+1]-d[l]) / (2.0*e[l]);
      r = pythag(g,1.0);
      g = d[m]-d[l]+e[l] / (g+sign(r,g));
      s = c = 1.0;
      p = 0.0;
      for(i=m-1;i>=l;i--)
        {
        f = s*e[i];
        b = c*e[i];
        e[i+1] = (r = pythag(f,g));
        if(r == 0.0)
          {
          d[i+1] -= p;
          e[m] = 0.0;
          break;
          }
        s = f/r;
        c = g/r;
        g = d[i+1]-p;
        r = (d[i]-g)*s+2.0*c*b;
        d[i+1] = g+(p=s*r);
        g = c*r-b;
        for(k=1;k<=dim;k++)
          {
          f=mat[k][i+1];
          mat[k][i+1] = s*mat[k][i]+c*f;
          mat[k][i] = c*mat[k][i]-s*f;
          }
        }
      if(r == 0.0 && i)
        continue;
      d[l] -= p;
      e[l] = g;
      e[m] = 0.0;
      }
    }
    while(m != l);
  }

} 

void householder(double **mat , int dim , double d[] , double e[])

/* Reduction de Householder Cf. Numerical Recipes in C */

{

double scale,hh,h,g,f;

for(i=dim;i>=2;i--)
  {
  l=i-1;
  h=scale=0.0;
  if(l > 1)
    {
    for(k=1;k<=l;k++)
      scale += fabs(mat[i][k]);
    if(scale == 0.0)
      e[i] = mat[i][l];
    else
      {
      for(k=1;k<=l;k++)
        {
        mat[i][k] /= scale;
        h += mat[i][k] * mat[i][k];
        }
      f = mat[i][l];
      g = (f > 0.0 ? -sqrt(h) : sqrt(h));
      e[i] = scale * g;
      h -= f*g;
      mat[i][l] = f-g;
      f = 0.0;
      for(j=1;j<=l;j++)
        {
        mat[j][i] = mat[i][j]/h;
        g = 0.0;
        for(k=1;k<=j;k++)
          g += mat[j][k] * mat[i][k];
        for(k=j+1;k<=l;k++)
          g += mat[k][j] * mat[i][k];
        e[j] = g/h;
        f += e[j] * mat[i][j];
        }
      hh = f / (h+h);
      for(j=1;j<=l;j++)
        {
        f = mat[i][j];
        e[j] = g = e[j] - hh * f;
        for(k=1;k<=j;k++)
          mat[j][k] -= (f*e[k]+g*mat[i][k]);
        }
      }
    }
  else
    e[i] = mat[i][l];
  d[i] = h;
  }

d[1] = 0.0;
e[1] = 0.0;
      
for(i=1;i<=dim;i++)
  {
  l=i-1;
  if(d[i])
    {
    for(j=1;j<=l;j++)
      {
      g = 0.0;
      for(k=1;k<=l;k++)
        g += mat[i][k] * mat[k][j];
      for(k=1;k<=l;k++)
        mat[k][j] -= g * mat[k][i];
      }
    }
  d[i] = mat[i][i];
  mat[i][i] = 1.0;
  for(j=1;j<=l;j++)
    mat[j][i] = mat[i][j] = 0.0; 
  }

}

long **matrix_l(int nrow, int ncol)

/* Allocation d'une matrice Cf. Numerical Recipes in C */

{

int ind1, ind2;

long **m;

m = (long **) malloc((size_t)((nrow+1)*sizeof(long*)));
if(!m)
  {
  printf("\nallocation failure 1 in matrix()");
  exit(0);
  }

m[1] = (long *) malloc((size_t)((nrow*(ncol+1))*sizeof(long)));
if(!m[1])
  {
  printf("\nallocation failure 2 in matrix()");
  exit(0);
  }

for(ind1=2;ind1<=nrow;ind1++)
  m[ind1]=m[ind1-1]+ncol;

return m;

}

double **matrix(int nrow, int ncol)

/* Allocation d'une matrice Cf. Numerical Recipes in C */

{

int ind1, ind2;

double **m;

m = (double **) malloc((size_t)((nrow+1)*sizeof(double*)));
if(!m)
  {
  printf("\nallocation failure 1 in matrix()");
  exit(0);
  }

m[1] = (double *) malloc((size_t)((nrow*(ncol+1))*sizeof(double)));
if(!m[1])
  {
  printf("\nallocation failure 2 in matrix()");
  exit(0);
  }

for(ind1=2;ind1<=nrow;ind1++)
  m[ind1]=m[ind1-1]+ncol;

return m;

}

long compute_cholesky(double **mat , int dim , double **vect)

/* Decomposition de Cholesky d'une matrice symetrique definie positive */

{

double sum;

for(i=1;i<=dim;i++)
  {
  for(j=i;j<=dim;j++)
    {
    for(sum=mat[i][j],k=i-1;k>=1;k--)
      sum -= mat[i][k] * mat[j][k];
    if(i == j)
      {
      if(sum <= 0.0)
        return(-1);
      vect[i][1] = sqrt(sum);
      }
    else
      mat[j][i] = sum / vect[i][1];
    }
  }

return(0);

}

void solve_sys(double **mat, int dim, double **vect, double **b , double **v)

/* Resolution d'un systeme lineaire */

{

double sum;

for(i=1;i<=dim;i++)
  {
  for(sum=b[i][1],k=i-1;k>=1;k--)
    sum -= mat[i][k] * v[k][1];
  v[i][1] = sum / vect[i][1];
  }

for(i=dim;i>=1;i--)
  {
  for(sum=v[i][1],k=i+1;k<=dim;k++)
    sum -= mat[k][i] * v[k][1]; 
  v[i][1] = sum / vect[i][1];
  }

}

void trans_mat(double **mat1 , double **mat2 , int dim1 , int dim2)
 
/* Calcul de la transposee d'une matrice */
 
{
 
for(i=1;i<=dim1;i++)
  for(j=1;j<=dim2;j++)
    mat2[j][i] = mat1[i][j];
 
}

double compute_trace(double *pmat , int dim)

/* Calcul de la trace d'une matrice */

{

int i;

double trace , *mat;

mat = pmat;
trace = *mat;

for(i=1;i<dim;i++)
  {
  mat += dim + 1;
  trace += *mat;
  }

return trace;

}

void add_mat(double **mat1 , double **mat2 , double **mat3 ,
                int dim1 , int dim2)
 
/* Addition de deux matrices */
 
{
 
int ind1, ind2;

for(ind1=1;ind1<=dim1;ind1++)
  for(ind2=1;ind2<=dim2;ind2++)
    mat3[ind1][ind2] = mat1[ind1][ind2] + mat2[ind1][ind2];

}
 
void sub_mat(double **mat1 , double **mat2 , double **mat3 ,
                int dim1 , int dim2)
 
/* Soustraction de deux matrices */
 
{
 
for(i=1;i<=dim1;i++)
  for(j=1;j<=dim2;j++)
    mat3[i][j] = mat1[i][j] - mat2[i][j];
 
}

void sub_mat1(double **mat1 , double **mat2 , int dim1 , int dim2)
 
/* Soustraction de deux matrices */
 
{
 
for(i=1;i<=dim1;i++)
  for(j=1;j<=dim2;j++)
    mat1[i][j] -= mat2[i][j];
 
}

void add_mat1(double **mat1 , double **mat2 , int dim1 , int dim2)
 
/* Addition de deux matrices, resultat stocke dans la premiere */
 
{
 
long ind1, ind2;

for(ind1=1;ind1<=dim1;ind1++)
  for(ind2=1;ind2<=dim2;ind2++)
    mat1[ind1][ind2] += mat2[ind1][ind2];
 
}

void scal_mat(double coeff , double **mat1 , double **mat2 ,
                 int dim1 , int dim2)
 
/* Multiplication d'une matrice par un scalaire */
 
{
 
for(i=1;i<=dim1;i++)
  for(j=1;j<=dim2;j++)
    mat2[i][j] = coeff * mat1[i][j];
 
}

void scal_mat1(double coeff , double **mat , int dim1 , int dim2)
 
/* Multiplication d'une matrice par un scalaire */
 
{
 
for(i=1;i<=dim1;i++)
  for(j=1;j<=dim2;j++)
    mat[i][j] *= coeff;
 
}

double gaussian(double *vect1 , double *vect2 , double dim)

/* Calcul d'une fonction noyau gaussienne */

{

long ind1;

double result=0.0, terme_partiel=0.0;

for(ind1=1; ind1<=dim; ind1++)
  {
  terme_partiel = (vect1[ind1] - vect2[ind1]);
  result -= (terme_partiel * terme_partiel);
  }

result /= (10.0 * (double)dim);

result = exp(result);

return result;

}

double prod_scal(double *vect1 , double *vect2 , long dim)

/* Produit scalaire de deux vecteurs (de meme dimension) */

{

long ind1;

double produit = 0.0;

for(ind1=1;ind1<=dim;ind1++)
  produit += vect1[ind1] * vect2[ind1];

return produit;
  
}

double polynomial(double *vect1 , double *vect2 , long dim)

/* Calcul d'une fonction noyau polynomiale */

{

long ind1;

double result=0.0;

for(ind1=1;ind1<=dim;ind1++)
  result += vect1[ind1] * vect2[ind1];
result += 1.0;
result = result * result;

return result;

}

double ker(int nat, double *vect1, double *vect2, long dim)

/* Calcul d'une fonction noyau */

{

double resultat;

switch(nat)
  {
  case 1: resultat = prod_scal(vect1, vect2, dim); break;
  case 2: resultat = gaussian(vect1, vect2, dim); break;
  case 3: resultat = polynomial(vect1, vect2, dim); break;
  default: printf("\n\nNoyau: %d inconnu...\n\n", nat); exit(0);
  }

return resultat;

}

void mult_mat(double **mat1 , double **mat2 , double **mat3 ,
                 int dim1 , int dim2 , int dim3)
 
/* Produit de deux matrices de tailles quelconques (compatibles) */
 
{
 
int ind1, ind2, ind3;

for(ind1=1;ind1<=dim1;ind1++)
  for(ind2=1;ind2<=dim3;ind2++)
    {
    mat3[ind1][ind2] = 0.0;
    for(ind3=1;ind3<=dim2;ind3++)
      mat3[ind1][ind2] += mat1[ind1][ind3] * mat2[ind3][ind2];
    }
 
}

void display_mat_l(long **mat , long dim1 , long dim2)

{

long ind1, ind2;

printf("\n");

for(ind1=1;ind1<=dim1;ind1++)
  {
  for(ind2=1;ind2<=dim2;ind2++)
    printf("%6d ", mat[ind1][ind2]);
  printf("\n");
  }

printf("\n");

}

void display_mat(double **mat , long dim1 , long dim2 , long par1 , long par2)

{

long ind1, ind2;

printf("\n");

for(ind1=1;ind1<=dim1;ind1++)
  {
  for(ind2=1;ind2<=dim2;ind2++)
    printf("%*.*f ",par1,par2,mat[ind1][ind2]);
  printf("\n");
  }

printf("\n");

}

void display_mat1(double **mat1 , double **mat2 , int dim1 , int dim2 ,
                  int par1 , int par2)
 
{
 
for(i=1;i<=dim1;i++)
  {
  for(j=1;j<=dim2;j++)
    printf("%*.*f %*.*f",par1,par2,mat1[i][j],par1,par2,mat2[i][j]);
  printf("\n");
  }
 
printf("\n");
 
}

void project(double **mat , int dim1 , int dim2)
 
/* Operateur de projection sur un cone */
 
{
 
for(i=1;i<=dim1;i++)
  for(j=1;j<=dim2;j++)
    if(mat[i][j] < 0.0)
      mat[i][j] = 0.0;
 
}

void copy(double **mat1 , double **mat2 , long dim1 , long dim2)

{

long id_ligne, id_col;

for(id_ligne=1;id_ligne<=dim1;id_ligne++)
  for(id_col=1;id_col<=dim2;id_col++)
    mat2[id_ligne][id_col] = mat1[id_ligne][id_col];

}

int compare(double **mat1 , double **mat2 , int dim1 , int dim2 , double seuil)

{

int id = true, end_mat = false;

i = j = 1;

while((id == true) && (end_mat == false))
  {
  if(fabs(mat1[i][j] - mat2[i][j]) > seuil)
    id = false;
  else
    {
    if(j < dim2)
      j++;
    else
      {
      if(i < dim1)
        {
        i++;
        j = 1;
        }
      else
        end_mat = true;
      }
    }
  }

return id;

}

void transpose_lignes(double **mat, int dim, int ligne1, int ligne2)

/* Transposition de deux lignes d'une matrice */

{

int indice;
double **tampon;

tampon = matrix(2, dim+1);

for(indice=1; indice<=dim; indice++)
  tampon[1][indice] = mat[ligne1][indice];

for(indice=1; indice<=dim; indice++)
  mat[ligne1][indice] = mat[ligne2][indice];

for(indice=1; indice<=dim; indice++)
  mat[ligne2][indice] = tampon[1][indice];

free(tampon);

}

double compute_det(double **mat, int dim)

/* Calcul du determinant d'une matrice */

{

long id1, id2, id3;
double **copie, pivot, quotient, negligeable = 1e-8, det;

copie = matrix(dim+1, dim+1);

for(id1=1; id1<=dim; id1++)
  for(id2=1; id2<=dim; id2++)
    copie[id1][id2] = mat[id1][id2];

det = 1.0;

for(id1=1; id1<dim; id1++)
  {
  pivot = copie[id1][id1];
  if(fabs(pivot) < negligeable)
    {
    id2 = id1;
    while((fabs(pivot) < negligeable) && (id2<dim))
      {
      id2++;
      pivot = copie[id2][id1];
      }
    if(fabs(pivot) >= negligeable)
      {
      transpose_lignes(copie, dim, id1, id2);
      det *= -1.0;
      }
    else
      {
      det = 0.0;
      free_matrix(copie);
      return det;
      }
    }

  for(id2=id1+1; id2<=dim; id2++)
    {
    quotient = copie[id2][id1] / copie[id1][id1];
    for(id3=id1; id3<=dim; id3++)
      copie[id2][id3] -= quotient * copie[id1][id3];
    }
  }
    
for(id1=1; id1<=dim; id1++)
  det *= copie[id1][id1];

free_matrix(copie);

return det;

}

int verif_inv(double **mat, double **mat_1, int dim)
 
/* Verification du calcul de l'inverse d'une matrice */
 
{
 
int id, id1, id2;
double **produit, **identite, negligeable = 1e-6;
 
produit = matrix(dim+1, dim+1);
identite = matrix(dim+1, dim+1);
 
for(id1=1; id1<= dim; id1++)
  for(id2=1; id2<= dim; id2++)
    identite[id1][id2] = (id1 == id2)? 1.0 : 0.0;
 
mult_mat(mat, mat_1, produit, dim, dim, dim);
 
id = compare(produit, identite, dim, dim, negligeable);
 
free(produit);
free(identite);
 
return(id);
 
}

void compute_inv(double **mat, double **mat_1, int dim)

/* Calcul de l'inverse d'une matrice carree quelconque */

{

long id1, id2, id3, id4;
double determinant, **mineur, negligeable = 1e-8;

mineur = matrix(dim, dim);

determinant = compute_det(mat, dim);
if(fabs(determinant) < pow(negligeable, (double) dim))
  {
  printf("\nLa matrice n'est pas inversible\n\n");
  return;
  }

for(id1=1; id1<=dim; id1++)
  for(id2=1; id2<=dim; id2++)
    {
    for(id3=1; id3<dim; id3++)
      for(id4=1; id4<dim; id4++)
        {
        if(id3 >= id1)
          {
          if(id4 >= id2)
            mineur[id3][id4] = mat[id3+1][id4+1];
          else
            mineur[id3][id4] = mat[id3+1][id4];
          }
        else
          {
          if(id4 >= id2)
            mineur[id3][id4] = mat[id3][id4+1];
          else
            mineur[id3][id4] = mat[id3][id4];
          }
        }
    mat_1[id2][id1] = compute_det(mineur, dim-1) / determinant;
    if((id1+id2)%2 == 1)
      mat_1[id2][id1] *= -1.0;
    }

if(verif_inv(mat, mat_1, dim) == false)
  {
  printf("\nEchec de l'inversion\n\n");
  display_mat(mat, dim, dim, 9, 5);
  display_mat(mat_1, dim, dim, 9, 5); 
  exit(0);
  }

free_matrix(mineur);

}