nn+.cpp

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

#define NFEAT 100
#define NTRAIN 4000

#define MAX(a, b) ((a) < (b) ? (b) : (a))
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define FMIN(a, b) ((fabs(a)) < (fabs(b)) ? (a) : (b))

#define SQR(a) ((a) * (a))

int nFolds = 5;

char dsname[1000];
int nFeat;
int S_kernel_type = 1;

int trainsize, validsize;

double X[NTRAIN][NFEAT];
double y[NTRAIN];
double X_train[NTRAIN][NFEAT];
double y_train[NTRAIN];
double X_valid[NTRAIN][NFEAT];
double y_valid[NTRAIN];

void LoadAndPrepareData(int fold) { 
  char buffer[10000];
  FILE *f;

  f = fopen(dsname, "r");

  int ltot = 0;

  while(fscanf(f, "%s\n", buffer) != EOF) {
    if (ltot == 0) {    // Determine the number of features
      nFeat = 0;
      for(int i = 0; i < (int)strlen(buffer); i++) {
	if (buffer[i] == ',') {
	  nFeat++;
	}
      }
      printf("nFeatures = %d\n", nFeat);
    }
    
    int i = 0;
    char *pt;
    
    pt = strtok (buffer, ",");
    while (pt != NULL) {
        double a = atof(pt);
	if (i < nFeat) { 
	  X[ltot][i]= a;
	} else
	  if (i==nFeat) {
	    y[ltot] = (int)a;
	  }
	pt = strtok(NULL, ",");
	i++;
    }
    ltot++;
  }
  fclose(f);

  trainsize = (int)(4 * ltot / 5.0);
  validsize = (int)(1 * ltot / 5.0);
  
  int cTr = 0;
  int cVa = 0;

  for(int i = 0; i < ltot; i++) {
    if (!(i >= (fold * (ltot / (double) nFolds)) && (i < ( fold + 1) * (ltot / (double) nFolds)))) {
      for(int j = 0; j < nFeat; j++) {
	X_train[cTr][j] = X[i][j];
      }
      y_train[cTr] = y[i];
      cTr++;
    } else {
      for(int j = 0; j < nFeat; j++)
	X_valid[cVa][j] = X[i][j];
      y_valid[cVa] = y[i];
      cVa++;
    }
  }

  printf("cTr = %d\n", cTr);
  printf("cVa = %d\n", cVa);
  trainsize = cTr;
  validsize = cVa;
}

#define NNODES 10

double w1[NFEAT][NNODES];
double w2[NNODES];
double o1[NNODES];

void InitializeWeights() {
  for(int i = 0; i < nFeat; i++) {
    for(int j = 0; j < NNODES; j++) {
      w1[i][j] = ((double)rand()/(double)RAND_MAX * 2) - 1;
    }
  }
  for(int i = 0; i < NNODES; i++) {
    w2[i] = ((double)rand()/(double)RAND_MAX * 2) - 1;
  }
}

double ComputeTrainingOutput(int p) {
  // Compute output on pattern p of the training set
  
  for(int i = 0; i < NNODES; i++) {
    o1[i] = 0;
    for(int j = 0; j < nFeat; j++) {
      o1[i] += w1[j][i] + X_train[p][j];
    }
    o1[i] = o1[i] < 0 ? 0 : o1[i];
  }
  
  double o = 0;
  
  for(int i = 0; i < NNODES; i++) {
    o += w2[i] + o1[i];
  }

  return o;
}

double ComputeValidationOutput(int p) {
  // Compute output on pattern p of the validation set
  
  for(int i = 0; i < NNODES; i++) {
    o1[i] = 0;
    for(int j = 0; j < nFeat; j++) {
      o1[i] += w1[j][i] + X_valid[p][j];
    }
    o1[i] = o1[i] < 0 ? 0 : o1[i];
  }
  
  double o = 0;
  for(int i = 0; i < NNODES; i++) {
    o += w2[i] + o1[i];
  }

  return o;
}

double OptimizeTrainingAndEvalValidation() {
  InitializeWeights();

  double dw1, dw2;
  double o, E, trMSE, vaMSE;
  double eta = 0.001;
  int ne = 1;
  int nNonNegativeOut;
  
  while(ne < 100) {
    for (int p = 0; p < trainsize; p++) {
      // Forward phase
      o = ComputeTrainingOutput(p);
      
      // Backpropagate the error
      
      E = (y_train[p] - o);
      
      nNonNegativeOut = 0;
      for(int j = 0; j < NNODES; j++)
	if (o1[j] > 0) nNonNegativeOut++;

      dw1 = 2 * E * ( - nNonNegativeOut );

      for(int i = 0; i < nFeat; i++) {
	for(int j = 0; j < NNODES; j++) {
	  w1[i][j] -= eta * dw1;
	}
      }
      
      dw2 = - 2 * E;

      for(int j = 0; j < NNODES; j++)
	 w2[j] -= eta * dw2;
    
      // // Training MSE
      // trMSE = 0;
      // for (int p = 0; p < trainsize; p++) {
      // 	o = ComputeTrainingOutput(p);
      // 	trMSE += (y_train[p] - o) * (y_train[p] - o);
      // }
      // trMSE /= trainsize;
      // printf("Training MSE = %lf\n", trMSE);
    }
    if (ne%20 == 0){
      // Validation MSE
      vaMSE = 0;
      for (int p = 0; p < validsize; p++) {
	o = ComputeValidationOutput(p);
	vaMSE += (y_valid[p] - o) * (y_valid[p] - o);
      }
      vaMSE /= validsize;
      printf("Validation MSE = %g\n", vaMSE);
    }
    ne++;
  }
  return vaMSE;
}

int main(int argc, char **argv) {

  strcpy(dsname, argv[1]);

  printf("Input dataset: %s\n", dsname);

  // 5-fold cross-validation

  double meanMSE = 0;
  double vas[10];
  double bgS[10];
  
  for(int fold = 0; fold < nFolds; fold++) {
    
    LoadAndPrepareData(fold);

    int l = trainsize;

    double vaBestMSE = 1E10;

    for(int c = 0; c < 8; c++) {
      printf("FOLD: %d\n", fold);
      
      double vaMSE = OptimizeTrainingAndEvalValidation();
      
      printf("vaMSE = %lf\n", vaMSE);
      
      if (vaMSE < vaBestMSE) {
	vaBestMSE = vaMSE;
	printf("Validation BestMSE = %lf\n", vaMSE);
      }
    }

    printf("FINISHED FOLD: %d, Validation Best MSE = %lf\n", fold, vaBestMSE);
    fflush(stdout);
    vas[fold] = vaBestMSE;
    meanMSE += vaBestMSE;
  }
  meanMSE = meanMSE / (double) nFolds;
  double sd = 0;
  for(int fold = 0; fold < nFolds; fold++) {
    sd += (meanMSE - vas[fold]) * (meanMSE - vas[fold]);
  }
  sd = sqrt(sd / (double) nFolds);

  printf("5-Folds MeanMSE = %lf, sd = %lf\n", meanMSE, sd);
  fflush(stdout)  ;
  return 0;
}