/* =======================================================
	Neural Network Classes for the NeXT Computer
	Written by: Ralph Zazula
					University of Arizona - Fall 1991
					zazula@pri.com (NeXT Mail)
==========================================================*/
#import "Neuron.h"
#import "Random.h"
#import <math.h>

#define	Np	16					// the number of points in the matching problem
#define 	N	Np*(Np-1)/2		// the number of neurons required to represent Np

id			random;				// random number generator instance
id			nList;				// list of the neurons
double	gamma = 1.0;		// "gamma" factor
double	alpha	= 0.995;		// temperature decay rate
double	p[Np][2];			// array to store the positions the points
double	d[Np][Np];			// array to store the distances between the points
double	T = 1.0;			   // the current temperature

//===========================================================
// function to return the index (from 0->N) of the neuron
// that represents the connection from point i->j
//
int index(i,j)
int	i,j;
{
	int	k;
	
	if(i>j) {				// swap them if reversed
		k = j;
		j = i;
		i = k;
	}
	return (int)(((float)Np - 3.0/2.0)*(float)i - (float)(i*i)/2.0 + j - 1);
}

//==========================================================
// function to return the distance between the points i and j
//
double dist(i,j)
int	i,j;
{
	return	sqrt((p[i][0] - p[j][0])*(p[i][0] - p[j][0]) +
					  (p[i][1] - p[j][1])*(p[i][1] - p[j][1]));
}

//===========================================================
// function to return the current value of the energy function
//
double H()
{
	double	h=0.0, x;
	
	int	i,j;

	
	for(i=0; i<Np; i++) {
		x = 0.0;
		for(j=0; j<Np; j++)
			if(i!=j)
				x += [[nList objectAt:index(i,j)] lastOutput];
		h += (gamma/2.0)*(1 - x)*(1 - x);
	}

	for(j=1; j<Np; j++)
		for(i=0; i<j; i++)
			h += d[i][j]*[[nList objectAt:index(i,j)] lastOutput];
	
	return h;
}	

//==========================================================
// function to return whether or not this node has flipped
// given the current temperature and the effect of the 
// flip on the energy funcion.
//
BOOL flip(n)
int	n;
{
	double	H1,H2,dH;
	
	//
	// get current energy-function, flip the node, get new energy-function
	// calculate dH
	//
	H1 = H();
	if([[nList objectAt:n] lastOutput])
		[[nList objectAt:n] setOutput:0];
	else 
		[[nList objectAt:n] setOutput:1];
							
	H2 = H();
	dH = H2 - H1;
	
	// 
	// flip back if the probablility was too low to flip this node
	// keep it flipped otherwise
	//
	if([random percent] > 1/(1 + exp(dH/T)))
		if([[nList objectAt:n] lastOutput])
			[[nList objectAt:n] setOutput:0];
		else
			[[nList objectAt:n] setOutput:1];

}

//===========================================================
// function to return the current length of all connections
//
double L()
{
	int	i,j;
	double l = 0.0;
	
	for(j=1; j<Np; j++)
		for(i=0; i<j; i++) 
			l += d[i][j]*[[nList objectAt:index(i,j)] lastOutput];
			
	return l;
}

//===========================================================

void main()
{
	int	i,j,count,n;
	
	nList = [[List alloc] init];				// list of neurons
	random = [[Random alloc] init];			// get a random number generator

	[random setSeeds:5335:7777:32197];
	//
	// generate Np random points
	//
	printf("generating points\n");
	for(i=0; i<Np; i++) {
		p[i][0] = [random percent];
		p[i][1] = [random percent];
	}
	// 
	// fill-in the distance array
	//
	printf("calculating distances\n");
	for(j=1; j<Np; j++)
		for(i=0; i<j; i++)
			d[i][j] = 
			d[j][i] = 
			 sqrt((p[i][0] - p[j][0])*(p[i][0] - p[j][0]) +
					(p[i][1] - p[j][1])*(p[i][1] - p[j][1]));
	//
	// generate Neurons
	//
	printf("generating neurons\n");
	for(i=0; i<N; i++) {
		[nList addObject:[[[Neuron alloc] init] setType:Binary]];
		[[nList lastObject] setOutput:[random randMax:1]];
//		printf("%f\n",[[nList lastObject] lastOutput]);
	}
	//
	// make initial connections
	//
	for(j=1; j<N; j++)
		for(i=0; i<j; i++) {
			[[nList objectAt:i] connect:[nList objectAt:j] withWeight:-gamma];
			[[nList objectAt:j] connect:[nList objectAt:i] withWeight:-gamma];
//			[[nList lastObject] step];
		}

/*
	for(i=0; i<Np; i++)
		printf("%u %f %f\n",i,p[i][0],p[i][1]);
	
	for(i=0; i<Np; i++)
		for(j=i+1; j<Np; j++)
			printf("%u->%u: %e\n",i,j,d[i][j]);
*/

	// start the network
	printf("starting total connection length: %f\n",L());
	printf("T    	 L  	   Energy     #of connections\n");
	for(count=0; count<200*N; count++)
	{
		//
		// pick a random node to update
		//
		n = [random randMax:N];
		[[nList objectAt:n] step];
		flip(n);
		if(!(count % N))	T = alpha*T;
		if(!(count % 100)) {
			n = 0;
			for(i=0; i<N; i++)
				if([[nList objectAt:i] lastOutput]) n++;
//			printf("count: %u  length:  %f  Temp: %f  Connections: %u\n",
//				count,L(),T,n);
			printf("%f %f %f %u\n",T,L(),H(),n);
		}
	}
	printf("points\n-----\n");
	for(i=0; i<Np; i++)
		printf("%u %f %f\n",i,p[i][0],p[i][1]);

	printf("distances\n-----\n");
	for(i=0; i<Np; i++)
		for(j=i+1; j<Np; j++)
				printf("%u->%u: %e\n",i,j,d[i][j]);

	printf("connections\n-----\n");
	for(i=0; i<Np; i++)
		for(j=i+1; j<Np; j++)
			if([[nList objectAt:index(i,j)] lastOutput])
				printf("%u->%u: %e\n",i,j,d[i][j]);

}
