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

#define INPUT_DIM	3	/* 入力データの次元数 */
#define N		500	/* 入力データ数 */
#define SOM_DIM		2	/* コードベクトルの次元数 */
#define L		6	/* 各次元のコードベクトル数 */
#define MAX_T		100	/* 最大繰り返し学習回数 */

#define rnd()		(((double)rand()/RAND_MAX-0.5)*2.0)
#define square(x)	((x)*(x))

/* 入力データ */
typedef double InputData[INPUT_DIM];
InputData X[N];

/* コードブック */
typedef struct {
	InputData *M;
	int n;
} CodeBookRec, *CodeBook;

/* 距離の定義 */
double Distance(InputData x, InputData m)
{
	register int i;
	double d;

	d=0.0;
	for(i=0; i<INPUT_DIM; i++)
		d+=square(x[i]-m[i]);
	return(d);
}

/* 学習係数関数 */
#define Lambda(t)	(1.0/(t))

/* コードベクトル間の距離 */
double DistanceSOM(int winner, int k)
{
	register int i;
	double d;
	int delta;

	d=0.0;
	for(i=0; i<SOM_DIM; i++) {
		delta=winner%L-k%L;
		d+=square(delta);
		winner/=L;
		k/=L;
	}
	return(sqrt(d));
}

/* 近傍関数 */
#define Phi(d)		(d<=1)

/* コードブックの初期化 */
CodeBook Initialize()
{
	register int i, j;
	CodeBook B;

	if((B=(CodeBook)malloc(sizeof(CodeBookRec)))==NULL)
		return(NULL);
	B->n=(int)pow(L, SOM_DIM);
	if((B->M=(InputData *)malloc(sizeof(InputData)*B->n))==NULL)
		return(NULL);
	/* コードブックの初期値設定 */
	for(i=0; i<B->n; i++)
		for(j=0; j<INPUT_DIM; j++)
			B->M[i][j]=rnd();
	return(B);
}

/* 最も近いコードベクトル番号を返す */
int Select(InputData x, CodeBook B)
{
	register int i;
	int min=-1;
	double d, min_d;

	for(i=0; i<B->n; i++) {
		d=Distance(x, B->M[i]);
		if(min<0 || d<min_d) {
			min=i;
			min_d=d;
		}
	}
	return(min);
}

/* SOM の学習 */
void Learn(InputData x, InputData m, double alpha)
{
	register int j;

	for(j=0; j<INPUT_DIM; j++)
		m[j]+=alpha*(x[j]-m[j]);
}

/* Winner とデータとの平均距離 */
double Delta(InputData X[], CodeBook B)
{
	register int i;
	double delta;
	int winner;

	delta=0;
	for(i=0; i<N; i++) {
		winner=Select(X[i], B);
		delta+=Distance(X[i], B->M[winner]);
	}
	return(delta/N);
}

main()
{
	register int t, i, j, k;
	CodeBook B;
	int winner;
	int *count;
	double r, delta, alpha;

	/* 入力データの生成(球のデータ) */
	for(t=0; t<N; t++) {
		r=0.0;
		for(j=0; j<INPUT_DIM; j++) {
			X[t][j]=rnd();
			r+=square(X[t][j]);
		}
		r=sqrt(r);
		for(j=0; j<INPUT_DIM; j++) {
			X[t][j]/=r;
		}
	}
	if((B=Initialize())==NULL)
		fprintf(stderr, "Cannot allocate Code Book\n");
	printf("t=0 delta=%lf\n", Delta(X, B));
	for(t=1; t<=MAX_T; t++) {
		for(i=0; i<N; i++) {
			winner=Select(X[i], B);
			for(k=0; k<B->n; k++) {
				alpha=Lambda(t)*Phi(DistanceSOM(winner, k));
				Learn(X[i], B->M[k], alpha);
			}
		}
		printf("t=%d delta=%lf\n", t, Delta(X, B));
	}
	/* 参照された回数のカウント */
	if((count=(int *)calloc(sizeof(int), B->n))==NULL) {
		fprintf(stderr, "Not enough memory\n");
		exit(1);
	}
	for(i=0; i<N; i++) {
		winner=Select(X[i], B);
		count[winner]++;
	}
	for(i=0; i<B->n; i++) {
		printf("%d: ", i);
		for(k=0; k<INPUT_DIM; k++)
			printf("%lf ", B->M[i][k]);
		printf("%d\n", count[i]);
	}
}