ploylinear.cpp

Locally weighted polynomial regression LWPR is a popular instance based al gorithm for learning c

  (1) It returns (int *num_terms) the number of terms
      associated with n input dimensions (for example, in linear
      regression with two inputs there are three terms. Quadratic
      regression with two inputs has 6 terms etc.  It ALWAYS does
      this.
  (2) IF *terms_vec is NULL it only does (1). But if
      *terms_vec is non-NULL it computes the terms values
      from input_vec and stores them in terms_vec. */
void linear_terms_function(double *input_vec,
			   int num_inputs,
			   int *num_terms,double *terms_vec)
{
  *num_terms = num_inputs + 1;
  if ( terms_vec != NULL )
  {
    int i;
    terms_vec[0] = 1.0;
    for ( i = 0 ; i < num_inputs ; i++ )
      terms_vec[i+1] = input_vec[i];
  }
}


/* Given inputs_mx, terms_mx, outputs_vec, query_vec, 
   weight_function, and terms_function are all defined,
   and memory has been allocated for weights_vec, 
   query_term_vec and beta_vec, this function will fill in...
     weights_vec (size num_points)
     query_term_vec (size num_terms)
     beta_vec (size num_terms)
   Returns the prediction: beta dot query_term. */
double lwpr_predict(double **inputs_mx,
                    double **terms_mx,
                    double *outputs_vec,
                    double *query_vec,
                    double *weights_vec,
                    double *query_term_vec,
                    double *beta_vec,
                    double (*weight_function)(double distance_sqd,
                                              double kernel_width),
		    void (*terms_function)(double *input_vec,
					   int num_inputs,
					   int *num_terms,
					   double *terms_vec),
		    int num_points,
		    int num_inputs,
		    int num_terms,
		    double kernel_width)
{
  int foo;
  double **wxtx_mx,  /* weighted X^transpose * X matrix */
    *wxty_vec,       /* weighted X^transpose * Y vector */  
    output;          /* prediction */

     /* filling in weights_vec */
  compute_weights(inputs_mx,query_vec,kernel_width,weight_function,
		  num_inputs,num_points, weights_vec);

     /* filling in query_term_vec */
  terms_function(query_vec, num_inputs, &foo, query_term_vec);
  
     /* allocating and computing wxtx_mx and wxty_vec */
  wxtx_mx = alloc_matrix(num_terms, num_terms);
  wxty_vec = alloc_vector(num_terms);
  compute_regression_matrices(terms_mx, outputs_vec, weights_vec, 
			      num_points, num_terms, wxtx_mx, wxty_vec);

     /* computes prediction and the regression 
	coefficients vecter beta_vec. */
  compute_beta(wxtx_mx, wxty_vec, num_terms, beta_vec);
  output = compute_predict(beta_vec, num_terms, query_term_vec);

    /* deallocates memory for wxtx_mx and wxty_vec and returns output */
  free_matrix(wxtx_mx, num_terms);
  free_vector(wxty_vec);
  return output;
}

/*  It allocates memory for and reads in in_mx and out_vec 
    from the filename, assuming that each row in the file 
    contains exactly 'number of input dimensions + 1' numbers, 
    separated by spaces  (it terminates with an error message 
    if not). It also returns the number of input dimensions
    and the number of data points. It is important not to 
    leave any blank lines at the beginning of the file. */
void load_data_from_filename(char *filename,  
			     double ***in_mx, 
			     double **out_vec,
			     int *num_inputs,
			     int *num_points)
{
  double **input_mx, *output_vec, foo;
  FILE *s;
  char c;
  int i, j;
  int prev; /* 0 indicates white-space; 1 indicates number character */
  int dims; /* number of input dimensions */
  int rows; /* number of data points */

  s = fopen(filename, "r");   /* attempts to open filename for reading */
  if(s==NULL)
  {
    printf("Error in openning %s.\n", filename);
    exit(-1);
  }
  
  dims = 0;    
  prev = 0;
  while((c=fgetc(s))!= '\n')     /* counts the number of dimensions */
    if(isspace(c) && prev == 1)
    {
      prev = 0;
      dims++;
    }
    else if(!isspace(c) && prev == 0)
      prev = 1;

  dims--;
  if(prev == 1)
    dims++;
 
  rows = 1;              /* counts the number of data points */
  while(fscanf(s, "%lf", &foo) != EOF) 
  {
    rows++;
    for(j=0; j<dims; j++)
      if(fscanf(s, "%lf", &foo) == EOF)
      {
	printf("Unexpected end of file while counting number of data points in %s.\n", filename);
	exit(-1);
      }
  }

  fclose(s);  /* closes filename */
 
  *num_inputs = dims;  /* sets reference variables */
  *num_points = rows;

  input_mx = alloc_matrix(rows, dims);  /* allocates memory */
  output_vec = alloc_vector(rows);

  s = fopen(filename, "r");  /* attempts to reopen filename for reading */
  if(s==NULL)
  {
    printf("Error in openning %s.\n", filename);
    exit(-1);
  }

  for(i=0; i<rows; i++)  /* inputs data from filename */
  { 
    for(j=0; j<dims; j++) 
      if(fscanf(s, "%lf", &(input_mx[i][j])) == EOF)
      {
	printf("Unexpected end of file while reading %s.\n", filename);
	exit(-1);
      }
    if(fscanf(s, "%lf", &(output_vec[i])) == EOF)
    {
      printf("Unexpected end of file while reading %s.\n", filename);
      exit(-1);
    }
  }
  fclose(s);   /* closes filename */

  *in_mx = input_mx;   /* sets reference variables */
  *out_vec = output_vec;
}


/* accepts the following command line arguments

   lwpr <datafile> <testfile> <kwidth>

  It reads in inputs_mx and outputs_vec from the datafile,
  assuming that each row in the datafile contains exactly the
  number of input dimensions + 1 numbers, separated by spaces 
  (it terminates with an error message if not). It is important
  not to leave any blank lines at the beginning of the files. 
  It similarly loads test_inputs_mx and test_outputs_vec from 
  testfile. It then does locally weighted linear regression 
  predictions on each row of test_inputs_mx. It prints to 
  standard output the predicted value and the true value from
  test_outputs_vec. At the end it prints the RMS prediction
  error.

  It frees up all its memory when done.
*/
void main(int argc,char *argv[])
{
  double **inputs_mx,  /* stores input matrix from datafile */
    **test_inputs_mx,  /* stores input matrix from testfile */
    **terms_mx,        /* terms used in regression corresponding to inputs_mx */
    *outputs_vec,      /* stores output vector from datafile */
    *test_outputs_vec, /* stores output vector from testfile */
    *weights_vec,      /* stores the weights for the current query point */
    *query_term_vec,   /* regression terms corresponding to the query point */
    *beta_vec,         /* coefficients of the regression */
    kernel_width,      /* user supplied kernel width */
    error,             /* root mean square error for predictions of testfile */
    val;               /* predicted output at the current query point */

  int num_inputs,      /* number of input dimensions */
    num_points,        /* number of points in datafile */
    test_num_inputs,   /* number of input dimensions in testfile; used for error checking */
    test_num_points,   /* number of points in testfile */
    num_terms,         /* number of columns in terms_mx */
    i;                 /* loop variable */
  
  if(argc != 4)   /* making sure there are enough arguments */
  {
    printf("Command line:\n<datafile> <testfile> <kwidth>\n\n");
    printf("where datafile contains the set of training data points,\n");
    printf("testfile contains the set of test data points,\n");
    printf("and kernel_width is the kernel width used in the locally weighted linear regression.\n");
    exit(-1);
  }

  kernel_width = atof(argv[3]);  /* extracts kernel width */
  
                  /* loading inputs_mx and test_inputs_mx from file names */ 
  load_data_from_filename(argv[1], &inputs_mx, &outputs_vec, 
			  &num_inputs, &num_points);
  load_data_from_filename(argv[2], &test_inputs_mx, &test_outputs_vec, 
			  &test_num_inputs, &test_num_points);

  printf("Number of input dimensions = %d.\n", num_inputs);
  printf("Number of training data points = %d.\n", num_points);
  printf("Number of test data points = %d.\n\n", test_num_points);
 
     /* checking that the number of input dimensions 
	equal in datafile and testfile */
  if(num_inputs != test_num_inputs)  
  {
    printf("Error: the number of input dimensions in your test and training files is not the same.\n");
    exit(-1);
  }

     /* allocating memory for weights_vec, terms_mx, query_term_vec,
	and beta_vec, and initializing num_terms */
  weights_vec = alloc_vector(num_points);
  terms_mx = compute_terms_mx(inputs_mx,num_inputs,num_points,
			      &num_terms, linear_terms_function);
  query_term_vec = alloc_vector(num_terms);
  beta_vec = alloc_vector(num_terms);

     /* loops through the test set, printing the predicted and correct
        output at each point and computes the root mean square error  
        of the predictions  */
  printf("Predicted Output------Actual Output\n");
  error = 0;
  for(i=0; i<test_num_points; i++)
  {
    val = lwpr_predict(inputs_mx,terms_mx,outputs_vec,
		       test_inputs_mx[i],weights_vec,
		       query_term_vec,beta_vec,gaussian_weight_function,
		       linear_terms_function, num_points, 
		       num_inputs, num_terms, kernel_width);
    printf("%14g %14g\n", val, test_outputs_vec[i]);
    error += SQR(val-test_outputs_vec[i]);
  }
  
  error /= test_num_points;
  error = sqrt(error);
  printf("\nRoot mean square error = %g.\n", error);
 

     /* freeing all allocated memory */
  free_matrix(inputs_mx, num_points);
  free_matrix(test_inputs_mx, test_num_points);
  free_matrix(terms_mx, num_points);
  free_vector(outputs_vec);
  free_vector(test_outputs_vec);
  free_vector(weights_vec);
  free_vector(query_term_vec);
  free_vector(beta_vec);
}