/*-------------------------------------------------------------

Copyright (C) 2000 Peter Clote. 
All Rights Reserved.

Permission to use, copy, modify, and distribute this
software and its documentation for NON-COMMERCIAL purposes
and without fee is hereby granted provided that this
copyright notice appears in all copies.


THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR
WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS
AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED
BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING
THIS SOFTWARE OR ITS DERIVATIVES.

-------------------------------------------------------------*/


/**********************************************

	Program: hopfield.c

This program implements Hopfield's discrete net
for the Traveling Salesman Problem, where N random cities
have real coordinates (x,y) in the unit square.
The approach is outlined in the text "Neural Networks
in Computer Intelligence" by LiMin Fu, McGraw Hill (1994),
especially pp 41--45. Updates are made at all nodes
simultaneously (ie synchronous updates).

	Programmed by: Peter Clote 
	Modified by: Christa Mauer
        Remodified by P. Clote spring 1998
       
**********************************************/



#include <stdio.h>
#include <math.h>       /* sqrt() prototype */
#include <stdlib.h>     /* def of RAND_MAX, rand() */
void error(char*);      /* prototype here so all functions can use */

#define STEPS 5     /* upper bound on number of steps to convergence */


#define A 500
#define B 500
#define C 50
#define D 500
#define N 5	/* N cities */
#define THETA 0.5



#define L(x)  ( (x) / (N)  )
#define R(x)  ( (x) % (N)  )
    /* For 0 <= x < N*N, x encodes its "left" and "right" parts i,j
       satisfying 0 <= i,j < N */

#define EPSILON 0.01    /* termination condition for main loop */


typedef double Array[N][N];
typedef struct 
    {
    double f;   /* first coordinate */
    double s;   /* second coordinate */
    } Point;        

typedef Point Pointarray[N];

/* copy(uold, unew) copies the NxN array unew to uold */
void copy ( Array uold, Array unew )
{
    int x,i;

    for (x=0;x<N;x++)
        for (i=0;i<N;i++)
            uold[x][i] = unew[x][i];
}

/*  successive returns 1 if x,y successive; i.e. either
y = (x+1) mod N or x = (y+1) mod N   */
int successive( int x, int y )
{
   if (  (y == ((x+1) % N)) ||  (x == ((y+1) % N)) )
      return 1;
   else
      return 0;
}

/* getpoints(points) initializes the Pointarray points of size N 
   with random points (f,s), where 0<=f,s<=1
 */
void getpoints ( Pointarray points )
{
    int i;

    for (i=0;i<N;i++)
        {
        points[i].f = (double) rand()/RAND_MAX;
        points[i].s = (double) rand()/RAND_MAX;
        }
}

/* getDistanceArray(points,dist) computes the distance array dist for the
   points in the point array points.
   for all x,y: dist[x][y] = euclidian_distance(points[x],points[y])
 */
void getDistanceArray ( Pointarray points, Array dist )
{
    double first, second;
    int x,y;
    
    for (x=0;x<N;x++)
      for (y=0;y<N;y++)
	{
	  first = points[x].f-points[y].f;
	  second = points[x].s-points[y].s;
	  dist[x][y] = sqrt(first*first + second*second);
	}
}

/* initializeOld(u) sets all entries in the NxN array u to 0.
 */
void initializeold ( Array u )
{
    int x,i;

    for (x=0;x<N;x++)
        for (i=0;i<N;i++)
            u[x][i] = 0;
}

/* initializenew(u) sets all entries in the NxN array u
  to random 0,1 values
*/
void initializenew ( Array u )
{
    int x,i;

    for (x=0;x<N;x++)
      for (i=0;i<N;i++)
	  u[x][i] = rand() % 2; 
}

/* difference(uold,unew) returns the maximal absolute difference
   between two corresponding entries in uold and unew

   returns max_0<=x,i<N |uold[x][i]-unew[x][i]|
*/
double difference ( Array uold, Array unew )
{
    int x,i;
    double temp,diff=0;

    for (x=0;x<N;x++)
        for (i=0;i<N;i++)
            {
            temp = abs(uold[x][i]-unew[x][i]);
            if ( temp > diff ) diff = temp;
            }
    return(diff);
}

double weight ( int a, int b, Array dist )
{
    /* weights of Hopfield net are taken as double */
    int x,y,i,j;

    x = L(a); i = R(a); y = L(b); j = R(b);
    	/* a = <x,i>, b = <y,j>  */
    if ( x != y )
        {
        if ( successive(i,j) == 1 )
            return ( -C - D * dist[x][y] );
        else if ( i != j )
            return ( -C );
        else        /* i == j */
            return ( -B-C );
        }
    else if ( i != j )  /* x == y */
        return ( -A-C );
    else
        return 0;
    /* if a == b, then weight(a,b,d) not called */
}


void update ( Array newo, Array dist )
{
    Array newMat;                /* Help */
    int i, j, n, nsquare;
    double weight(int a, int b, Array dist);
    int sum;
    
    nsquare = N*N;
    for (i=0;i<nsquare;i++)
        {
        newMat[L(i)][R(i)] = 0;                 
        }
    for (i=0;i<nsquare;i++)
        {
        sum = 0;
        for (j=0;j<nsquare;j++) 
            if ( i != j ) sum += weight(j,i,dist)*newo[L(j)][R(j)];
        if ( sum > THETA )
            newMat[L(i)][R(i)] = 1;               
                        /* Temporary storage of new values */
        else if ( sum < THETA )
            newMat[L(i)][R(i)] = 0;
        else 
            newMat[L(i)][R(i)] = newo[L(i)][R(i)];
        }
    for (i=0;i<nsquare;i++)
        {
        newo[L(i)][R(i)] = newMat[L(i)][R(i)];
        }
}

main ()
{
    /********** function prototypes ***********/
    void initializenew ( Array );
    void initializeold ( Array );
    void display ( Array );
    void getpoints ( Pointarray points );
    void getDistanceArray ( Pointarray, Array );
    double difference ( Array, Array);
    void update ( Array, Array );
    void copy ( Array, Array );
    double weight ( int, int, Array );
    void printpoints ( Pointarray );

    /*********** variables *******************/
    Array uold, unew, dist;
    Pointarray points;
    double diff=1, deltaT = 0.01; /* initial difference, time interval */
    int i;
     

    /*********** initialize stuff *******************/
    getpoints(points); /* initialize points at random */
    getDistanceArray(points,dist); /* compute distance matrix */
    initializeold(uold);    /* uold is 0 */
    initializenew(unew);    /* unew is random */
    
    printf("Initialization:\n");
    printf("Points\n");
    printpoints(points);
    printf("Uold\n");
    display(uold);
    printf("Unew\n");
    display(unew);
    printf("\n");

    /**************** process ************************/

    diff = difference(uold,unew);
    for (i=0;i<STEPS && diff>=EPSILON;i++)
        {
        copy(uold,unew);
        update(unew,dist);
        diff = difference(uold,unew);
        printf("Step:%d \t Difference %f\n",i,diff); 
        printf("Uold\n");
	display(uold);
        printf("Unew\n");
        display(unew);
	printf("\n");
        }
    if ( diff >= EPSILON )
        printf("No convergence after %d steps\n",STEPS);
    else
        {
        printf("Convergence after %d steps\n",i-1);
        printf("Difference is %lf\n",diff);
        printf("Uold\n");
        display(uold);
        printf("Unew\n");
        display(unew);        
        }


}



/************************ routines for debugging ***********/



void printpoints ( Pointarray points )
{
    int i;

    for (i=0;i<N;i++)
        printf("%lf %lf\n",points[i].f,points[i].s);
}


void display ( Array u )
{
    int i,x;
    for (x=0;x<N;x++)
        {
        for (i=0;i<N;i++)
			printf("%8.5lf ",u[x][i]);
        printf("\n");
        }
}


void error ( char *s )
{
    fprintf(stderr,"%s\n",s);
    exit(1);
}