/*

   Original Author: Shawn P. Legrand      15-Sep-88

	This program simulates a simple neural network by  calculating the 
    inner product of a vector and a matrix where the matrix T is symmetric
    {T(i,j) = T(j,i)}, dilute (more 0's than 1's), and T(i,i)=0 (to prevent
    oscillations and chaotism). Reports are sent to a report file showing the
    input and output matrices and the energy level of the matrix after each
    iteration.	

	This program is based on information contained in a book by Ed Rietman,
    'Experiments in Artificial Neural Networks', published by TAB.

    Copyright 1988 	Shawn P. Legrand, CCP 		All Rights Reserved
	
*/

#include <stdio.h>

/*
 * Random number generator -
 * adapted from the FORTRAN version 
 * in "Software Manual for the Elementary Functions"
 * by W.J. Cody, Jr and William Waite.
*/

static long int iy = 100001;

sran(seed)
long seed;
{
        iy = seed;
}

double ran()
{
        iy *= 125;
        iy -= (iy/2796203) * 2796203;
        return (double) iy/ 2796203.0;
}

/*
 * Main program logic.
*/

main()
{
    int neurons, io, info, vector, matrix, i, j, iterate, sigma=0;
    long seed;
    float energy=0;
    double r;
    int t[100][100], u[100], v[100];
    FILE *fp;

/* Open new CLI device */

    if ((fp = fopen("Neuron.output","w")) == NULL) {
       printf("Error opening Neuron.output!!!\n");
       exit(10);
    }   

/* Retrive initial values */

    printf("Input the random seed: ");
    scanf("%ld",&seed);
    sran(seed);
     
    printf("Enter the number of neurons (100 maximum): ");
    scanf("%d",&neurons);						       	
    printf("Input the threshold value (0 to 2 are reasonable values): ");
    scanf("%d",&io);
    printf("Enter the value of the information (0 to 1 is a good value): ");
    scanf("%d",&info);
    printf("Do you want to enter the input vector yourself (1/Yes, 0/No)? ");
    scanf("%d",&vector);
    printf("Do you want to enter the T matrix (1/Yes, 0/No)? ");
    scanf("%d",&matrix);

/* Initialize the weight matrix (matrix T) */    

    if (!matrix) {
    	for (i=0; i<neurons; i++) {
	    for (j=0; j<neurons; j++) {
	        if (i==j) {
		    t[i][j] = 0;
		} else {
                    r=ran();
		    if (r<0.8) {
		        t[i][j] = 0;
		    } else {
		        t[i][j] = 1;
		    }
		}    
	    }
	}
    } else {
        for (i=0; i<neurons; i++) {
	    for (j=0; j<neurons; j++) {
	        printf("T(%d,%d): ",i,j);
		scanf("%d",&t[i][j]);
	    }
        }
    }

/* Output initial values */
    
    fprintf(fp,"Seed = %ld\n",seed); 
    fprintf(fp,"Threshold = %d\n",io);
    fprintf(fp,"Information = %d\n",info);
    fprintf(fp,"T matrix:\n");
    for (i=0; i<neurons; i++) {
        for (j=0; j<neurons; j++) {
	    fprintf(fp,"%d",t[i][j]);
	}
	fprintf(fp,"\n");
    }
    fprintf(fp,"\n\n");    
    
/* Retrieve input vector (vector u) */

    if (!vector) {
        for (i=0; i<neurons; i++) {
	    r=ran();
	    if (r<0.5) {
	        u[i]=0;
	    } else {
	        u[i]=1;
	    }
	}
    } else {
        for (i=0; i<neurons; i++) {
	    fprintf(fp,"Input u(%d): ",i);
	    fscanf(fp,"%d",&u[i]);
	}
    }

/* Perform calculations, Report input and output vector, Energy level, 
   Iterate using feedback from prior calculation */    

    for (iterate=0; iterate<8; iterate++) {
        for (i=0; i<neurons; i++) {
	    for (j=0; j<neurons; j++) {
	        sigma+=t[i][j]*u[j];
	    }
	    sigma+=info;
	    if (sigma>io) {
	        sigma=1;
	    } else {
	        sigma=0;
	    }
	    v[i]=sigma;
	    sigma=0;
	}
	fprintf(fp,"Iteration %d:\n",iterate);
	fprintf(fp,"Input vector U:\n");
	for (i=0; i<neurons; i++) {
	    fprintf(fp,"%d",u[i]);
	}
	fprintf(fp,"\nOutput vector V:\n");
	for (i=0; i<neurons; i++) {
	    fprintf(fp,"%d",v[i]);
	}
	fprintf(fp,"\n");
	fprintf(fp,"            Energy: ");
	for (i=0; i<neurons; i++) {
	    energy+=(float) (u[i]*v[i]);
	}
	energy*=-0.5;
	fprintf(fp,"%f\n\n",energy);
	for (i=0; i<neurons; i++) {
	    u[i]=v[i];
	}
	energy=0;
    }
    
}