/*
 * rbf.C
 *
 * For Problem 3 in Tutorial 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.5; x[0][1] = 1.2;
  x[1][0] = 0.0; x[1][1] = -0.4; 
  x[2][0] = 0.5; x[2][1] = 0.7;
  x[3][0] = 1.2; x[3][1] = 0.5; 
  x[4][0] = -0.1; x[4][1] = 1.5;
  x[5][0] = 0.3; x[5][1] = 0.8; 
  x[6][0] = 1.0; x[6][1] = -0.5; 
  x[7][0] = 0.4; x[7][1] = -1.0;
  x[8][0] = 0.8; x[8][1] = -0.2;
  x[9][0] = 0.2; x[9][1] = -0.1;

// assign indexe of cluster to each input pattern

  cluster[0] = 0; cluster[1] = 0; cluster[2] = 0;
  cluster[3] = 1; cluster[4] = 1; cluster[5] = 1;
  cluster[6] = 2; cluster[7] = 2; cluster[8] = 2; cluster[9] = 2;


// find cluster centers

  for (i=0; i<3; i++) { 
    count[i] = 0; 
    for (j=0; j<2; j++) 
      sum[i][j] = 0.0;
  }
  
  for (i=0; i<10; i++) {
    color = cluster[i];
    count[color]++;
    for (j=0; j<2; j++)
      sum[color][j] += x[i][j];
  }
	
  for (i=0; i<3; i++) 
    for (j=0; j<2; j++)
      mean[i][j] = sum[i][j]/count[i];
  
  fprintf(stdout, "mean values \n");

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

  // find covariance matrices 
  
  for (k=0; k<3; k++)
    for (i=0; i<2; i++)
      for (j=0; j<2; j++) 
	cov[k][i][j] = 0.0;

  for (k=0; k<10; k++) {
    color = cluster[k];
    for (i=0; i<2; i++)
      for (j=0; j<2; j++) 
	cov[color][i][j] += (x[k][i] - mean[color][i])*(x[k][j] - mean[color][j]);
  }

  for (k=0; k<3; k++)
    for (i=0; i<2; i++)
      for (j=0; j<2; j++)
	cov[k][i][j] /= (count[k] - 1.0);
  
  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<10; 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<10; 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[3][0]*0.1 + r[6][0]*0.8 + r[8][0]*0.5;
  double r1d = r[0][1]*(-0.5) + r[3][1]*0.1 + r[6][1]*0.8 + r[8][1]*0.5;
  double r2d = r[0][2]*(-0.5) + r[3][2]*0.1 + r[6][2]*0.8 + r[8][2]*0.5;
  
  double r0r1 = r[0][0]*r[0][1] + r[3][0]*r[3][1] + r[6][0]*r[6][1] + r[8][0]*r[8][1];
  double r1r2 = r[0][2]*r[0][1] + r[3][2]*r[3][1] + r[6][2]*r[6][1] + r[8][2]*r[8][1];
  double r0r2 = r[0][0]*r[0][2] + r[3][0]*r[3][2] + r[6][0]*r[6][2] + r[8][0]*r[8][2];

  double r0r0 = SQR(r[0][0]) + SQR(r[3][0]) + SQR(r[6][0]) + SQR(r[8][0]);
  double r1r1 = SQR(r[0][1]) + SQR(r[3][1]) + SQR(r[6][1]) + SQR(r[8][1]);
  double r2r2 = SQR(r[0][2]) + SQR(r[3][2]) + SQR(r[6][2]) + SQR(r[8][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]);
  
}