/*
 * rbf.C
 *
 * For Example 3 in Section 9
 *
 */

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

#define SQR(x) ((x)*(x))

main()
{
  int i, j, k;
  
  double det;
  double x[10][3], mean[3][2], sum[3][2], cov[3][2][2], inv[3][2][2];
  double r[10][3], d[10];
  int cluster[10], count[3], color;
  int flips;

 
// initialize the seven patterns

  x[0][0] = -1.0; x[0][1] = 0.2;
  x[1][0] = 0.2; x[1][1] = 1.0; 
  x[2][0] = 0.5; x[2][1] = -1.2;
  x[3][0] = 1.2; x[3][1] = -0.2; 
  x[4][0] = 0.8; x[4][1] = 2.1;

	
   mean[0][0]=-1.0; mean[0][1]=0.0;
   mean[1][0]=0.5; mean[1][1]=0.4;
   mean[2][0]=0.2; mean[2][1]=-1.5;
  

  // find covariance matrices 
  
   cov[0][0][0] = 1.0;
   cov[0][0][1] = 0.4;
   cov[0][1][0] = 0.5;
   cov[0][1][1] = 1.2;

   cov[1][0][0] = 2.0;
   cov[1][0][1] = 0.1;
   cov[1][1][0] = 0.2;
   cov[1][1][1] = 1.2;

   cov[2][0][0] = 0.75;
   cov[2][0][1] = 0.3;
   cov[2][1][0] = 0.4;
   cov[2][1][1] = 1.2;
  
  fprintf(stdout,"covariance matrices \n");

  for (i=0; i<2; i++) {
    for (k=0; k<3; k++) 
      fprintf(stdout, "|%6.2f %6.2f |   ", cov[k][i][0], cov[k][i][1]);
    fprintf(stdout, "\n");
  }
  fprintf(stdout, "\n");

  // find inverse of covariance matrices
  for (k=0; k<3; k++) {
    det = cov[k][0][0]*cov[k][1][1] - cov[k][0][1]*cov[k][1][0];
    inv[k][0][0] = cov[k][1][1]/det;
    inv[k][1][1] = cov[k][0][0]/det;
    inv[k][0][1] = -cov[k][0][1]/det;
    inv[k][1][0] = -cov[k][1][0]/det;
  }
  
  fprintf(stdout, "inverse of covariance \n");

  for (i=0; i<2; i++) {
    for (k=0; k<3; k++) 
      fprintf(stdout, "|%6.2f %6.2f |   ", inv[k][i][0], inv[k][i][1]);
    fprintf(stdout, "\n");
  }
  fprintf(stdout, "\n");
  
  // find the hidden layer outputs for the input patterns
  for (k=0; k<5; k++) {
    r[k][0] = exp(-0.5*(SQR(x[k][0]-mean[0][0])*inv[0][0][0] +
			2.0*(x[k][0]-mean[0][0])*(x[k][1]-mean[0][1])*inv[0][0][1] +
			SQR(x[k][1]-mean[0][1])*inv[0][1][1]));
    r[k][1] = exp(-0.5*(SQR(x[k][0]-mean[1][0])*inv[1][0][0] +
			2.0*(x[k][0]-mean[1][0])*(x[k][1]-mean[1][1])*inv[1][0][1] +
			SQR(x[k][1]-mean[1][1])*inv[1][1][1]));
    r[k][2] = exp(-0.5*(SQR(x[k][0]-mean[2][0])*inv[2][0][0] +
			2.0*(x[k][0]-mean[2][0])*(x[k][1]-mean[2][1])*inv[2][0][1] +
			SQR(x[k][1]-mean[2][1])*inv[2][1][1]));
  }

  fprintf(stdout,"hidden layer activities for inputs \n");

  for (k=0; k<5; k++)
    fprintf(stdout, " (%6.2f, %6.2f %6.2f )\n ", r[k][0], r[k][1], r[k][3]);
  fprintf(stdout,"\n");



  // find summations for LMS algorithm

  double r0d = r[0][0]*(0.5) + r[1][0]*(-1.0) + r[2][0]*1.2 + r[3][0]*(-0.75) + r[4][0]*1.0;
  double r1d = r[0][1]*(0.5) + r[1][1]*(-1.0) + r[2][1]*1.2 + r[3][1]*(-0.75) + r[4][1]*1.0;
  double r2d = r[0][2]*(0.5) + r[1][2]*(-1.0) + r[2][2]*1.2 + r[3][1]*(-0.75) + r[4][2]*1.0;
  
  double r0r1 = r[0][0]*r[0][1] + r[1][0]*r[1][1] + r[2][0]*r[2][1] + r[3][0]*r[3][1]+ r[4][0]*r[4][1];
  double r1r2 = r[0][2]*r[0][1] + r[1][2]*r[1][1] + r[2][2]*r[2][1] + r[3][2]*r[3][1] + r[4][2]*r[4][1];
  double r0r2 = r[0][0]*r[0][2] + r[1][0]*r[1][2] + r[2][0]*r[2][2] + r[3][0]*r[3][2] + r[4][0]*r[4][2];

  double r0r0 = SQR(r[0][0]) + SQR(r[1][0]) + SQR(r[2][0]) + SQR(r[3][0]) + SQR(r[4][0]);
  double r1r1 = SQR(r[0][1]) + SQR(r[1][1]) + SQR(r[2][1]) + SQR(r[3][1]) + SQR(r[4][1]);
  double r2r2 = SQR(r[0][2]) + SQR(r[1][2]) + SQR(r[2][2]) + SQR(r[3][2])+ SQR(r[4][2]);

  fprintf(stdout, "summations \n");
  fprintf(stdout, "r0d  =%6.2f,   r1d  =%6.2f,   r2d  =%6.2f \n", r0d, r1d, r2d);
  fprintf(stdout, "r0r1 =%6.2f,   r1r2 =%6.2f,   r0r2 =%6.2f \n", r0r1, r1r2, r0r2);
  fprintf(stdout, "r0r0 =%6.2f,   r1r1 =%6.2f,   r2r2 =%6.2f \n\n", r0r0, r1r1, r2r2);


  // NOW SOLVE LINEAR EQUATIONS USING MATRICES
  
  double A[3][3]; 

  A[0][0] = r0r0; A[0][1] = r0r1; A[0][2] = r0r2;
  A[1][0] = r0r1; A[1][1] = r1r1; A[1][2] = r1r2;
  A[2][0] = r0r2; A[2][1] = r1r2; A[2][2] = r2r2;

// find the Determinent

  det = A[0][0]*(A[1][1]*A[2][2] - A[1][2]*A[2][1]);
  det -= A[0][1]*(A[1][0]*A[2][2] - A[1][2]*A[2][0]);
  det += A[0][2]*(A[1][0]*A[2][1] - A[1][1]*A[2][0]);

// alert for singular matrices
  double eps = 0.00001;
  if (det >= 0.0 && det < eps) det = eps;
  else if (det < 0.0 && det > -eps) det = -eps;
    
// find the inverse of A;
  double invA[3][3];
  
  invA[0][0] = (A[1][1]*A[2][2] - A[2][1]*A[1][2])/det;
  invA[0][1] = (A[2][1]*A[0][2] - A[0][1]*A[2][2])/det;
  invA[0][2] = (A[0][1]*A[1][2] - A[1][1]*A[0][2])/det;
  invA[1][0] = (A[1][2]*A[2][0] - A[1][0]*A[2][0])/det;
  invA[1][1] = (A[0][0]*A[2][2] - A[2][0]*A[0][2])/det;
  invA[1][2] = (A[1][0]*A[0][2] - A[0][0]*A[1][2])/det;
  invA[2][0] = (A[1][0]*A[2][1] - A[1][1]*A[2][0])/det;
  invA[2][1] = (A[0][1]*A[2][0] - A[0][0]*A[2][1])/det;
  invA[2][2] = (A[0][0]*A[1][1] - A[0][1]*A[1][0])/det;

  // find weights
  double w[3];
  for (i = 0; i < 3; i++) 
    w[i] = invA[i][0]*r0d + invA[i][1]*r1d + invA[i][2]*r2d;

  fprintf(stdout, "weights \n");
  fprintf(stdout, "%6.2f %6.2f %6.2f \n", w[0], w[1], w[2]);
  
}