xtal_v51.c

统计模式识别算法包

/*****************************************************************************/
/*    XTAL - A Multimethod Program for Multivariate Regression               */
/*                                                                           */
/*           Don Gehring - April 14,1994 - Version 5.1                       */
/*****************************************************************************/



/* ------------------------------------------------------------------------ */
/* These parameters can be changed to alternate settings                    */
#define SPLINE 1   /* spline parameter for MARS */
#define RETRIES 9  /* simulated annealling retries for ANN */
#define ITERATIONS 50  /* iteration parameter for CTM */
char ANNINITIALIZATION[] = "ANNEAL NOREGRESS"; /*  ANN init procedure */
/* ------------------------------------------------------------------------ */

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>

double pardub ( char **str ) ;
int cntcol( char **str ) ;
void star_bar_plot(FILE *pfile, double y , double yscale );
int lookup_cmd(char *cmd);
int one_param(int *a);
int two_param( int *a, int *b);
void ann1();
void ctm1();
void mars1();
void ppr1();
void knn1();

char parameter_line[241];
char *para_ptr;
char trainname[16],testname[16] ;
int trn_samp_num,tst_samp_num,predx ;
char out_file_name[] = "rms.out";
char maxerr_name[] = "max.out" ;
double response_mean_rms,response_mean_variance;
double tst_diff,tst_mean ;
double trn_min,trn_max,tst_min,tst_max ;
double trnx_min[11],tstx_min[11],trnx_max[11],tstx_max[11];
double response_min,response_max,scale ;
double xmin[11],xmax[11],scalex[11];
char temptrainname[] = "tempor.trn";
char temptestname[] = "tempor.tst";
char tem2trainname[] = "tempor2.trn";
char tem2testname[] = "tempor2.tst" ;


int main()
{
  char filenamecontrol[] = "filename.con" ;
  char operationcontrol[] = "operate.con" ;
  int i,j,k,n;
  char *str ;
  char **pstr ;
  double mse_sum,rms,y,ytemp,y_est,diff,x ;
  int flag ,command_num;
  char line[241] , command[5];
  FILE *trainfile , *testfile , *outfile ,*maxerr;
  FILE *namefile , *operfile;
  double response_sum,response_mean,response;
  FILE *temptrainfile, *temptestfile;
  FILE *tem2trainfile, *tem2testfile;
  char s[30] ;
  /* create file to write program output into
     this file is first opened and reclosed after use to 
     allow the user to examine it during program execution.
     Its contents allows the user to determine how far the 
     program is from completion*/
  outfile = fopen(out_file_name, "w");
  fclose(outfile);
  /* create file to write the worst case error values into
     this file is first opened and reclosed to erase it */
  maxerr = fopen(maxerr_name, "w");
  fclose(maxerr);
  /* open filename control file */
  namefile = fopen(filenamecontrol,"r+");
  if(namefile == NULL) {
    printf("\nCan't open %s - job terminated\n",filenamecontrol);
    return 0 ;
  }

  /* open operation control file */
  operfile = fopen(operationcontrol,"r+");
  if(operfile == NULL) {
    printf("\nCan't open %s - job terminated\n",operationcontrol);
    return 0 ;
  }

  /******* MAIN OUTER LOOP BEGINS HERE ****************/
  do{
    /* open file to write program output into
       this file is first opened and reclosed after use to 
       allow the user to examine it during program execution.
       Its contents allows the user to determine how far the 
       program is to completion */  
    outfile = fopen(out_file_name,"a");
  loop1:  /* jump back to here if illegal input line */
    
    /* read a line from the filename control file */
    flag = fscanf(namefile,"%s %s",trainname,testname);
    if(flag == EOF) break;  /* no more data in file - break do loop */

    /* check to make sure there were two filenames read */
    if(flag != 2) {
      fprintf(outfile,"Wrong number of filenames on a line in filename control file");
      goto loop1;
    }
    trainfile = fopen(trainname,"r");
    if(trainfile == NULL) {
      fprintf(outfile,"Can't open file %s \n",trainname);
      goto loop1;
    }
    
    trn_samp_num = 0;
    response_sum = 0.0 ;
    trn_min = 9999999999.0 ;
    trn_max = -9999999999.0 ;
    for (i=1 ; i<=10 ; i++) {
      trnx_min[i] = 99999999999.9;
      tstx_min[i] = 99999999999.9;
      trnx_max[i] = -99999999999.9;
      tstx_max[i] = -99999999999.9;
    }

    do{
      flag = fscanf(trainfile,"%[^\n]\n",line);
      if(flag == EOF) 
	break;
      if(trn_samp_num == 0) {
	str = line;
	pstr = &str ;
	predx = cntcol(pstr) - 1 ;
      }
      str = line;
      pstr = &str ;
      for(i=1; i<=predx; i++) {
	x = pardub(pstr);
	if(x < trnx_min[i]) trnx_min[i] = x;
	if(x > trnx_max[i]) trnx_max[i] = x;
      }

      /* get the response variable from train file input line */
      /* use it to compute the response mean */
      response = pardub(pstr);
      response_sum += response ;
        
      if(response < trn_min) 
	trn_min = response;
      if(response > trn_max) 
	trn_max = response;
      trn_samp_num++;
    }while(1);
    fclose(trainfile);
    if((predx <= 0) || (predx >= 10)) {
      fprintf(outfile,"Wrong numbers of predictor variables in training file %s\n",trainfile);
      goto loop1 ; 
    }
   
    /* compute the response mean */
    response_mean = response_sum / (double)trn_samp_num ;
    testfile = fopen(testname,"r");
    if(testfile == NULL) {
      fprintf(outfile,"Can't open file %s \n",testname);
      goto loop1 ;
    }
    tst_samp_num = 0;
    mse_sum = 0.0;
    tst_mean = 0.0;
    tst_min = 9999999999.0 ;
    tst_max = -9999999999.0 ;

    do{
      flag = fscanf(testfile,"%[^\n]\n",line);
      if(flag == EOF) 
	break;
      tst_samp_num++;
      str = line;
      pstr = &str ;
      for(i=1 ; i<=predx ; i++) {
	x = pardub(pstr);
	if(x < tstx_min[i]) tstx_min[i] =x;                                                    
	if(x > tstx_max[i]) tstx_max[i] = x;
      }

      /* get the response variable from test file input line */
      y = pardub(pstr);
      diff = response_mean - y ;
      mse_sum += diff * diff ;
      tst_mean += y ;
      if(y < tst_min) tst_min = y ;
      if(y > tst_max) tst_max = y ;
    }while(1);

    fclose(testfile);
    response_mean_variance = mse_sum/(tst_samp_num+1) ;
    response_mean_rms = sqrt(response_mean_variance);
    tst_mean=tst_mean/(tst_samp_num);
    tst_diff = tst_max-tst_min ;

    /* find scaling factor for each predictor variable */
    for (i=1 ; i<=predx ; i++) {
      if(trnx_min[i]<tstx_min[i]) 
	xmin[i] = trnx_min[i];
      else 
	xmin[i] = tstx_min[i];
      if(trnx_max[i] > tstx_max[i]) 
	xmax[i] = trnx_max[i];
      else 
	xmax[i] = tstx_max[i];
      scalex[i] = xmax[i] - xmin[i] ;
      
    }
    /************* rescale training file so all response values are 
  between .1 and .9 *****************************************/
    
    if(trn_min < tst_min) 
      response_min = trn_min;
    else 
      response_min = tst_min;

    if(trn_max > tst_max) 
      response_max = trn_max;
    else 
      response_max = tst_max ;
    scale = response_max - response_min ;

    trainfile = fopen(trainname,"r");
    temptrainfile = fopen(temptrainname,"w");
    tem2trainfile = fopen(tem2trainname,"w");

    do{
      flag= fscanf(trainfile,"%[^\n]\n",line);
      if(flag==EOF) 
	break;
      str = line;
      pstr = &str;

      for( n=1 ; n<=predx ; n++) {
	x = -1.0 +( 2.0*(pardub(pstr) - xmin[n]))/scalex[n]; 
	fprintf(temptrainfile,"%14.12f  ",x);
	fprintf(tem2trainfile,"%14.12f  ",x);
      }
      
      ytemp = pardub(pstr);
      y = .1 +( (ytemp-response_min)*.8)/scale;
      fprintf(temptrainfile,"%14.12f\n",y);
      fprintf(tem2trainfile,"%14.12f\n",ytemp);
    }while(1);
      
    fclose(trainfile);
    fclose(temptrainfile);
    fclose(tem2trainfile);
    testfile = fopen(testname,"r");
    tem2testfile = fopen(tem2testname,"w");
    temptestfile = fopen(temptestname,"w");
    
    do{
      flag= fscanf(testfile,"%[^\n]\n",line);
      if(flag==EOF) break;
      str = line;
      pstr = &str;

      for( n=1 ; n<=predx ; n++) {
	x = -1.0 +( 2.0*(pardub(pstr) - xmin[n]))/scalex[n]; 
	fprintf(tem2testfile,"%14.12f  ",x);
	fprintf(temptestfile,"%14.12f  ",x);
      }

      ytemp = pardub(pstr);
      y = .1 +( (ytemp-response_min)*.8)/scale;
      fprintf(temptestfile,"%14.12f\n",y);
      fprintf(tem2testfile,"%14.12f\n",ytemp);
    }while(1);
    fclose(testfile);
    fclose(temptestfile);
    fclose(tem2testfile);
   
    for(i=1; i<=80 ; i++) 
      fprintf(outfile,"_");
    fprintf(outfile,"\n\nTraining file %s opened with %d samples and %2d predictor variables (mean=%f)\n",trainname,trn_samp_num,predx,response_mean);
    fprintf(outfile,"Test file %s opened with %d samples (mean=%f)\n",testname,tst_samp_num,tst_mean);
    fprintf(outfile,"Zeroth order model (y=%f) has RMS error=%f on the test set.\n", response_mean,response_mean_rms);
    for(i=1; i<=80 ; i++) 
      fprintf(outfile,"_");
    fprintf(outfile,"\n");    
    fclose(outfile);

    /* minmaxout is the output file that shows worst case error */
    maxerr = fopen(maxerr_name,"a");
    for(i=1; i<=80 ; i++) 
      fprintf(maxerr,"_");
    fprintf(maxerr,"\n\nTraining file %s has response MIN = %f\n",trainname,trn_min);
    fprintf(maxerr,"Training file %s has response MAX = %f\n",trainname,trn_max);
    fprintf(maxerr,"Test file %s has response MIN = %f\n",testname,tst_min);
    fprintf(maxerr,"Test file %s has response MAX = %f\n",testname,tst_max);
    fprintf(maxerr,"The test response MIN-MAX difference is %f\n",tst_diff);
    for(i=1; i<=80 ; i++) 
      fprintf(maxerr,"_");
    fprintf(maxerr,"\n");
    fclose(maxerr);

    /************ INNER CONTROL FILE LOOP BEGINS HERE ********************/
    fseek(operfile, 0 , 0 );   /* reset operation file ptr to zero */
    do{
    loop2: /* if error in control line return here */
      flag = fscanf(operfile,"%[^\n]\n",line);
      if(flag == EOF) 
	break; /* control file now empty */
      if(strlen(line) < 4) {
	fprintf(outfile,"Too few characters in command line\n");
	goto finish;
      }

      strncpy(command,line,4);
      command[4] = '\0' ;
      command_num = lookup_cmd(command);
      if(command_num == 0){
	fprintf(outfile,"Illegal command or syntax in operation file\n");
	goto finish;
      }

      strcpy(parameter_line, line+4) ;
      para_ptr = parameter_line ;
      switch (command_num){
      case 1 : 
	ann1(); 
	break;
      case 2 : 
	ctm1(); 
	break;
      case 3 : 
	mars1(); 
	break;
      case 4 : 
	ppr1(); 
	break;
      case 5 : 
	knn1(); 
	break;
      }
    }while(1);
  }while(1);
 finish:
  
  strcpy(s,"rm ");
  strcat(s,temptrainname);
  system(s);

  strcpy(s,"rm ");
  strcat(s,temptestname);
  system(s);

  strcpy(s,"rm ");
  strcat(s,tem2trainname);
  system(s);

  strcpy(s,"rm ");
  strcat(s,tem2testname);
  system(s);

  fclose(operfile);
  fclose(namefile);
  return 1 ;
}

/*_________________________________________________________________________*/
/*                                                                         */
/*  PARSDUBL - ParseDouble routine to parse a double from a string         */
/*                                                                         */
/*_________________________________________________________________________*/
double pardub ( char **str )
{
  double num ;
  while (! (isdigit ( **str )  
	    ||  (**str == '-')  
	    ||  (**str == '.') 
	    || (**str == '+') )) {
    if (**str)
      ++(*str) ;
    else
      return 0. ;
  }
  num = atof ( *str ) ;
  while ((isdigit ( **str )
	  ||(**str=='-')
	  ||(**str=='.')
	  ||(**str=='e')
	  ||(**str=='E') 
	  || (**str == '+') ))
    ++(*str) ;
  
  return num ;
}

/*_________________________________________________________________________*/
/*                                                                         */
/*  CNTCOL - routine to count the number of numbers in a line of chars     */
/*                                                                         */
/*_________________________________________________________________________*/
int cntcol ( char **str )
{
  int i = 0 ;
  do{
    while (! (isdigit ( **str )  
	      ||  (**str == '-')  
	      ||  (**str == '.') 
	      || (**str == '+') )) {
      if (**str)
	++(*str) ;
      else
	return i ;
    }
    i++;
    while (isdigit ( **str )
	   ||(**str=='-')
	   ||(**str=='.')
	   ||(**str=='e')
	   ||(**str=='E')
	   ||(**str=='+') )
      ++(*str) ;
  }while(1);
}

void star_bar_plot(FILE *pfile, double y , double yscale ) {
#define COLUMNS 80
  int jj,ypos ;
  int columns = COLUMNS ;
  char matrix[COLUMNS+1] ;
  for(jj=1 ; jj<=columns ; jj++){
    if( jj%5 != 0 ) matrix[jj] = '.' ;
    else matrix[jj] = '|' ;
  }

  ypos =1 + (int)( (double)(columns-1) * y  / yscale ) ;
  if(ypos > columns) ypos = columns;
  if(ypos < 1) ypos = 1 ;
  ypos = columns-ypos+1 ;

  for(jj=1 ; jj<=ypos ; jj++) 
    matrix[jj] = '*' ;

  for(jj=1 ; jj<=columns ; jj++) 
    fprintf(pfile,"%c",matrix[jj]);
  fprintf(pfile,"\n");
  return ;
}

int lookup_cmd(char *cmd){
  if(strcmp(cmd,"ANN1") == 0) 
    return 1 ;
  if(strcmp(cmd,"CTM1") == 0) 
    return 2 ;
  if(strcmp(cmd,"CTMR") == 0) 
    return 2 ;
  if(strcmp(cmd,"MRS1") == 0) 
    return 3 ;
  if(strcmp(cmd,"MARS") == 0) 
    return 3 ;
  if(strcmp(cmd,"PPR1") == 0) 
    return 4 ;
  if(strcmp(cmd,"PRPR") == 0) 
    return 4 ;
  if(strcmp(cmd,"KNN1") == 0) 
    return 5 ;
  if(strcmp(cmd,"K-NN") == 0) 
    return 5 ;
  return 0;
}

int one_param(int *a){
  int num;
  while( ! ( isdigit(*para_ptr)) ) {
    if(*para_ptr) 
      ++(para_ptr);
    else 
      return 0 ;
  }
  num = atoi( para_ptr );
  while ( isdigit(*para_ptr) ) 
    ++(para_ptr) ;

  *a = num ;
  return 1 ;
}

int two_param(int *a ,int *b){
  int num;