discrim.c

source codes for ltr

	  }
	}
	beta[virtual_class][att] = temp;
      }
    }
  }
  return beta;
}

static double *inner\
(DomainInfo *domain, double ***STATIS, double **beta, int nratt, int nrcla, double *class_freq)
{
  register int att, cl, virtual_class;
  double temp, *alfa;
 
  alfa = dvector(1, NR_CLASS);

  for(cl = 1; cl <= nrcla; cl++) {
    virtual_class = CLASS_USED[cl];
    if (virtual_class) {
      temp = 0.0;
      for(att = 1; att <= nratt; att++) {
	if (CiTypeAttr(domain, att) == continuous) {
	  temp += beta[virtual_class][att] * STATIS[att][cl][1];
	}
	else  {
	  temp += beta[virtual_class][att] * STATIS[att][cl][1+NValsAttr(domain, att)];
	}
      }
      alfa[virtual_class] = temp;
    }
  }

  for(cl = 1; cl <= nrcla; cl++) {
    virtual_class = CLASS_USED[cl];
    if (virtual_class) {
      alfa[virtual_class] =  HALF * alfa[virtual_class] + log((double) class_freq[cl] / class_freq[0]);
    }
  }
  for(cl = 1; cl <= nrcla; cl++)    alfa[cl] -= alfa[NR_CLASS];
  return alfa;
}

static void corrige_beta(double **beta, int nratt, int nrcla)
{
  register int i, j, cl;

  for(i = 1; i <= nrcla; i++) {
    cl = CLASS_USED[i];
    if (cl)
      for(j = 1; j <= nratt; j++)
	beta[cl][j] -= beta[NR_CLASS][j];
  }
}

/********************************************************/
/*	Singular Value Decomposition	A = U.W.Vt	*/
/*	Input 	a[1..m][1..n]				*/
/*	Output	a[1..m][1..n]	as U			*/
/*		w[1..n]					*/
/*		v[1..n][1..n]				*/
/********************************************************/
void svd_treshold(double tresh)
{
  TRESH = tresh;
}

int svdcmp(double **a, int m, int n, double *w, double **v)
{
  int flag, i, its, j, jj, k, l, nm;
  double anorm, c, f, g, h, s, scale, x, y, z, *rv1;

  rv1 = dvector(1, n);
  g = scale = anorm = 0.0;
  for(i=1;i<=n;i++) {
    l = i + 1;
    rv1[i]=scale * g;
    g = s = scale = 0.0;
    if (i <= m) {
      for(k=i;k<=m;k++)	scale += fabs(a[k][i]);
      if (scale) {
	for(k=i;k<=m;k++){
	  a[k][i] /= scale;
	  s += SQR(a[k][i]);
	}
	f = a[i][i];
	g = -SIGN(sqrt(s), f);
	h = f * g - s;
	a[i][i] = f- g;
	for(j=l;j<=n;j++) {
	  for(s=0.0, k= i;k<=m;k++)	s += a[k][i]*a[k][j];
	  f = s/h;
	  for(k=i;k<=m;k++)	a[k][j] += f*a[k][i];
	}
	for(k=i;k<=m;k++)	a[k][i] *= scale;
      }
    }
    w[i] = scale * g;
    g = s = scale = 0.0;
    if (i <= m && i != n) {
      for(k=l; k <=n ; k++)	scale += fabs(a[i][k]);
      if (scale) {
	for(k=l; k<=n; k++) {
	  a[i][k] /= scale;
	  s += SQR(a[i][k]);
	}
	f = a[i][l];
	g = -SIGN(sqrt(s), f);
	h = f * g - s;
	a[i][l] = f - g;
	for(k=l;k<=n;k++)	rv1[k]=a[i][k]/h;
	for(j=l;j<=m;j++) {
	  for(s=0.0, k=l; k<=n; k++)	s+=a[j][k]*a[i][k];
	  for(k=l;k<=n;k++)		a[j][k] += s*rv1[k];
	}
	for(k=l; k<=n; k++)	a[i][k] *= scale;
      }
    }
    anorm = MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
  }
  for(i=n;i>=1;i--) {
    if(i < n) {
      if (g) {
	for(j=l; j<=n; j++)	v[j][i]=(a[i][j]/a[i][l])/g;
	for(j=l; j<=n; j++) {
	  for(s=0.0,k=l;k<=n;k++)	s += a[i][k]*v[k][j];
	  for(k=l;k<=n;k++)	v[k][j] += s*v[k][i];
	}
      }
      for(j=l;j<=n;j++)	v[i][j] = v[j][i] = 0.0;
    }
    v[i][i] = 1.0;
    g = rv1[i];
    l = i;
  }
  for(i=MIN(m,n);i>=1;i--) {
    l = i+1;
    g=w[i];
    for(j=l; j<=n; j++)	a[i][j] = 0.0;
    if(g) {
      g = 1.0 / g;
      for(j=l; j<=n; j++) {
	for(s=0.0,k=l; k <= m; k++)	s += a[k][i]*a[k][j];
	f = (s/a[i][i])*g;
	for(k=i; k<=m; k++)	a[k][j] += f*a[k][i];
      }
      for(j=i;j<=m;j++)	a[j][i] *= g;
    }
    else	for(j=i;j<=m;j++)	a[j][i] = 0.0;
    ++a[i][i];
  }
  for(k=n;k>=1;k--) {
    for(its=1; its<=30; its++) {
      flag = 1;
      for(l=k;l>=1;l--) {
	nm = l-1;
	if((double)(fabs(rv1[l])+anorm) == anorm) {
	  flag = 0;
	  break;
	}
	if ((double)(fabs(w[nm])+anorm) == anorm)	break;
      }
      if (flag) {
	c = 0.0;
	s = 1.0;
	for(i=l;i<=k;i++) {
	  f = s * rv1[i];
	  rv1[i] = c * rv1[i];
	  if ((double)(fabs(f)+anorm) == anorm) break;
	  g = w[i];
	  h = pythag(f,g);
	  w[i] = h;
	  h = 1.0 / h;
	  c = g*h;
	  s = -f*h;
	  for(j=1; j<=m; j++) {
	    y=a[j][nm];
	    z=a[j][i];
	    a[j][nm] = y*c+z*s;
	    a[j][i] = z*c-y*s;
	  }
	}
      }
      z=w[k];
      if (l == k) {
	if (z < 0.0) {
	  w[k] = -z;
	  for(j=1;j<=n;j++)	v[j][k] = -v[j][k];
	}
	break;
      }
      if(its == 30) {
	fprintf(stderr,"E)SVD: No Convergence\n");
	return 0;
      }
      x = w[l];
      nm = k-1;
      y = w[nm];
      g = rv1[nm];
      h = rv1[k];
      f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
      g = pythag(f, 1.0);	
      f = ((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
      c = s = 1.0;
      for(j=l; j<=nm; j++) {
	i = j+1;
	g = rv1[i];
	y = w[i];
	h = s*g;
	g = c*g;
	z = pythag(f,h);
	rv1[j] = z;
	c = f/z;
	s = h/z;
	f = x*c+g*s;
	g = g*c-x*s;
	h = y*s;
	y *= c;
	for(jj=1; jj<=n; jj++) {
	  x = v[jj][j];
	  z = v[jj][i];
	  v[jj][j] = x*c+z*s;
	  v[jj][i] = z*c-x*s;
	}
	z = pythag(f, h);
	w[j] = z;
	if (z) {
	  z = 1.0/z;
	  c = f*z;
	  s = h*z;
	}
	f = c*g+s*y;
	x = c*y-s*g;
	for(jj = 1; jj<=m; jj++) {
	  y = a[jj][j];
	  z = a[jj][i];
	  a[jj][j] = y*c+z*s;
	  a[jj][i] = z*c-y*s;
	}
      }
      rv1[l] = 0.0;
      rv1[k] = f;
      w[k] = x;
    }
  }
  free_dvector(rv1, 1, n);
  return 1;
}

static void svbksb(double **u, double *w, double **v, int m, int n, double *b, double *x)
{
  int i, j, jj;
  double s, *tmp, wmin, wmax = 0.0;
  
  for(j = 1; j <= n; j++)	if(w[j] > wmax)	wmax = w[j];
  wmin = wmax * TRESH;
  for(j = 1; j <= n; j++)	if(w[j] < wmin)	w[j] = 0.0;
  
  tmp = dvector(1, n);
  for(j = 1; j <= n; j++) {
    s = 0.0;
    if (w[j]) {
      for(i = 1; i <= m; i++)	s += u[i][j] * b[i];
      s /= w[j];
    }
    tmp[j] = s;
  }
  for(j = 1; j <= n; j++) {
    s = 0.0;
    for(jj = 1; jj <= n; jj++)	s += v[j][jj] * tmp[jj];
    x[j] = s;
  }
  free_dvector(tmp, 1, n);
}
/********************************************************/
/*	Inversao de matrizes por SVD	      	        */
/*	Input:	a[1..m][1..n]				*/
/*	Output	y[1..m][1..n]				*/
/********************************************************/
static double **svd_inv(double **a, int m, int n)
{
  double **y = NULL, **v, *w;
  double *col = dvector(1, n);
  double wmax, wmin;
  int i, j;
  
  v = dmatrix(1, n, 1, n);
  w = dvector(1, n);

  if (svdcmp(a, m, n, w, v)) {
    wmax = w[1];
    for(j = 1; j <= n; j++) if (w[j] > wmax) wmax = w[j];
    wmin = wmax * TRESH;
    for(j = 1; j <= n; j++) if (w[j] < wmin) w[j] = 0.0; 
    y = dmatrix(1, m, 1, n);
    for(j = 1; j <= m; j++) {
      for(i = 1; i <= n; i++)	col[i] = 0.0;
      col[j] = 1.0;
      svbksb(a, w, v, m, n, col, col);
      for(i = 1; i <= n; i++)	y[i][j] = col[i];
    }
    free_dvector(w, 1, n);
    free_dmatrix(v, 1, n, 1, n);
  }
  free_dvector(col, 1, n);
  return(y);
}

static double pythag(double a, double b)
{
  double absa, absb;
  absa = fabs(a);
  absb = fabs(b);
  if (absa > absb)      return(absa*sqrt(1.0+SQR(absb/absa)));
  return(absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
}