/* 
 * Necker Cube alternation simulation, syncronous algorithm
 * Author: ASAKAWA Shinichi <asakawa@twcu.ac.jp>
 * $Id: Necker-alt.c,v 1.1 2002/04/29 07:56:16 asakawa Exp $
 *
 * Configuration of neurons:
 * 
          0 ------- 1              8 ------- 9
         /         /|	          /|         |
        2---------3 |	        10 |      11 |
        |         | |	         | |         |
        | 4       | 5	         |12---------13
        |         |/ 	         |/         / 
        6---------7  	        14---------15 
*/

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

#define N_neuron 16

int FROM=0;
int TO=10000;
int EPOCH=1;

double sigmoid_slope = -2.0;
double Inhibition = -2.0;
double Excitation =  1.0;
double Noise_range = 10.0;

/* 
 * Definition of a neuron
 * Ninp strands for a number of input from other neurons.
 */
struct NEURON {
    double output;
    double state;
    int fire_cnt;
    double thres;
    int    Ninp;
    double *inp_wgt;
    double (*f)();
    struct NEURON **inp_neuron;
};
typedef struct NEURON neuron;

neuron Necker[N_neuron];


double linear(double value)
{
	return value;
}

double sigmoid(double value)
{
	return 1.0 / ( 1.0 + exp(sigmoid_slope * value) );
}

double formal(double value)
{
    if ( value > 0.0 ){
	return 1.0;
    } else {
	return 0.0;
    }
}

/* return random number from 0 to 1 */
double drand(void)
{
	return ((double)(rand() % RAND_MAX) / (double)RAND_MAX);
}

/* return random number from -1 to 1 */
double drand2(void)
{
	return (drand() * 2.0 - 1.0);
}

void initialize_neuron(neuron *a, int N)
{
	int i;

	a->thres = 0.0;
	a->state = drand();
	a->fire_cnt = 0;
	a->Ninp = N;
	a->inp_wgt = (double *)malloc(sizeof(double) * (N+1));
	for (i=0; i<=N; i++) {
		a->inp_wgt[i] = 0.0;
	}
	a->inp_neuron = (neuron **)malloc(sizeof(neuron *) * (N+1));
	a->f = sigmoid;
	a->output = a->f(a->state);
}

void calc_neuron(neuron *a)
{
    double wrk = 0.0;
    int i;

    for(i=0; i< a->Ninp; i++){
	wrk += a->inp_wgt[i] * a->inp_neuron[i]->output;
    }
    wrk -= a->thres;
    a->state = wrk;
    a->output = a->f(a->state);
}

void Print_net(int iter)
{
	printf("### Iteration : %d\n", iter);
	printf("   %3.1f ------ %3.1f                 %3.1f ----- %3.1f\n",
	       Necker[0].output,Necker[1].output,Necker[8].output,Necker[9].output);
	printf("    /         /|                  /|         |\n");
        printf(" %3.1f ------ %3.1f|               %3.1f |     %3.1f |\n",
	       Necker[2].output,Necker[3].output,Necker[10].output,Necker[11].output);
        printf("   |         | |                 | |         |\n");
        printf("   | %3.1f     | %3.1f	         |%3.1f ----- %3.1f\n",
	       Necker[4].output,Necker[5].output,Necker[12].output,Necker[13].output);
        printf("   |         |/                  |/         / \n");
        printf(" %3.1f ------ %3.1f                %3.1f ------ %3.1f \n\n",
	       Necker[6].output,Necker[7].output,Necker[14].output,Necker[15].output);
}

int main_loop(void)
{
	int i, iter = 0;

	while ( iter <= TO ){
	    for (i=0;i<N_neuron;i++){
		if (Necker[i].output > 0.8) {
		    Necker[i].fire_cnt++;
		}
		if (Necker[i].fire_cnt > 10) {
		    Necker[i].state = -30.0;
		    Necker[i].output = 0.0;
		    Necker[i].fire_cnt = -10;
		} 
		if (Necker[i].fire_cnt < 0 ) {
		    Necker[i].state = -30.0;
		    Necker[i].output =  0.0;
		    Necker[i].fire_cnt++;
		} else {
		    Necker[i].state += Noise_range * drand();
		}
	    }
	    for (i=0;i<N_neuron;i++){
		calc_neuron(&Necker[i]);
	    }
	    if ( (iter % EPOCH) == 0 ) {
		Print_net(iter);
	    }
	    iter++;
	}
	return EXIT_SUCCESS;
}

void Make_Necker_Cube_Network(void)
{
	int i;

	for (i=0; i<16; i++) {
		initialize_neuron(Necker+i, 4);
	}

	/* Mutal inhibitary connections */
	for (i=0; i<8; i++) {
		Necker[i].inp_neuron[0] = &Necker[i+8];
		Necker[i].inp_wgt[0] = Inhibition;

		Necker[i+8].inp_neuron[0] = &Necker[i];
		Necker[i+8].inp_wgt[0] = Inhibition;
	}

	/* Excitatory connections */
	/* left side */
	Necker[0].inp_neuron[1] = &Necker[1]; Necker[0].inp_wgt[1] = Excitation;
	Necker[0].inp_neuron[2] = &Necker[2]; Necker[0].inp_wgt[2] = Excitation;
	Necker[0].inp_neuron[3] = &Necker[4]; Necker[0].inp_wgt[3] = Excitation;
				     
	Necker[1].inp_neuron[1] = &Necker[0]; Necker[1].inp_wgt[1] = Excitation;
	Necker[1].inp_neuron[2] = &Necker[3]; Necker[1].inp_wgt[2] = Excitation;
	Necker[1].inp_neuron[3] = &Necker[5]; Necker[1].inp_wgt[3] = Excitation;
				     
	Necker[2].inp_neuron[1] = &Necker[0]; Necker[2].inp_wgt[1] = Excitation;
	Necker[2].inp_neuron[2] = &Necker[3]; Necker[2].inp_wgt[2] = Excitation;
	Necker[2].inp_neuron[3] = &Necker[6]; Necker[2].inp_wgt[3] = Excitation;
				     
	Necker[3].inp_neuron[1] = &Necker[1]; Necker[3].inp_wgt[1] = Excitation;
	Necker[3].inp_neuron[2] = &Necker[2]; Necker[3].inp_wgt[2] = Excitation;
	Necker[3].inp_neuron[3] = &Necker[7]; Necker[3].inp_wgt[3] = Excitation;
	
	Necker[4].inp_neuron[1] = &Necker[0]; Necker[4].inp_wgt[1] = Excitation;
	Necker[4].inp_neuron[2] = &Necker[5]; Necker[4].inp_wgt[2] = Excitation;
	Necker[4].inp_neuron[3] = &Necker[6]; Necker[4].inp_wgt[3] = Excitation;
				     
	Necker[5].inp_neuron[1] = &Necker[1]; Necker[5].inp_wgt[1] = Excitation;
	Necker[5].inp_neuron[2] = &Necker[4]; Necker[5].inp_wgt[2] = Excitation;
	Necker[5].inp_neuron[3] = &Necker[7]; Necker[5].inp_wgt[3] = Excitation;
				     
	Necker[6].inp_neuron[1] = &Necker[2]; Necker[6].inp_wgt[1] = Excitation;
	Necker[6].inp_neuron[2] = &Necker[4]; Necker[6].inp_wgt[2] = Excitation;
	Necker[6].inp_neuron[3] = &Necker[7]; Necker[6].inp_wgt[3] = Excitation;
				     
	Necker[7].inp_neuron[1] = &Necker[3]; Necker[7].inp_wgt[1] = Excitation;
	Necker[7].inp_neuron[2] = &Necker[5]; Necker[7].inp_wgt[2] = Excitation;
	Necker[7].inp_neuron[3] = &Necker[6]; Necker[7].inp_wgt[3] = Excitation;

	/* right side */
	Necker[0+8].inp_neuron[1] = &Necker[1+8]; Necker[0+8].inp_wgt[1] = Excitation;
	Necker[0+8].inp_neuron[2] = &Necker[2+8]; Necker[0+8].inp_wgt[2] = Excitation;
	Necker[0+8].inp_neuron[3] = &Necker[4+8]; Necker[0+8].inp_wgt[3] = Excitation;
	         		       
	Necker[1+8].inp_neuron[1] = &Necker[0+8]; Necker[1+8].inp_wgt[1] = Excitation;
	Necker[1+8].inp_neuron[2] = &Necker[3+8]; Necker[1+8].inp_wgt[2] = Excitation;
	Necker[1+8].inp_neuron[3] = &Necker[5+8]; Necker[1+8].inp_wgt[3] = Excitation;
	         		       
	Necker[2+8].inp_neuron[1] = &Necker[0+8]; Necker[2+8].inp_wgt[1] = Excitation;
	Necker[2+8].inp_neuron[2] = &Necker[3+8]; Necker[2+8].inp_wgt[2] = Excitation;
	Necker[2+8].inp_neuron[3] = &Necker[6+8]; Necker[2+8].inp_wgt[3] = Excitation;
	         		       
	Necker[3+8].inp_neuron[1] = &Necker[1+8]; Necker[3+8].inp_wgt[1] = Excitation;
	Necker[3+8].inp_neuron[2] = &Necker[2+8]; Necker[3+8].inp_wgt[2] = Excitation;
	Necker[3+8].inp_neuron[3] = &Necker[7+8]; Necker[3+8].inp_wgt[3] = Excitation;

	Necker[4+8].inp_neuron[1] = &Necker[0+8]; Necker[4+8].inp_wgt[1] = Excitation;
	Necker[4+8].inp_neuron[2] = &Necker[5+8]; Necker[4+8].inp_wgt[2] = Excitation;
	Necker[4+8].inp_neuron[3] = &Necker[6+8]; Necker[4+8].inp_wgt[3] = Excitation;
	         		            
	Necker[5+8].inp_neuron[1] = &Necker[1+8]; Necker[5+8].inp_wgt[1] = Excitation;
	Necker[5+8].inp_neuron[2] = &Necker[4+8]; Necker[5+8].inp_wgt[2] = Excitation;
	Necker[5+8].inp_neuron[3] = &Necker[7+8]; Necker[5+8].inp_wgt[3] = Excitation;
	         		            
	Necker[6+8].inp_neuron[1] = &Necker[2+8]; Necker[6+8].inp_wgt[1] = Excitation;
	Necker[6+8].inp_neuron[2] = &Necker[4+8]; Necker[6+8].inp_wgt[2] = Excitation;
	Necker[6+8].inp_neuron[3] = &Necker[7+8]; Necker[6+8].inp_wgt[3] = Excitation;
	         		            
	Necker[7+8].inp_neuron[1] = &Necker[3+8]; Necker[7+8].inp_wgt[1] = Excitation;
	Necker[7+8].inp_neuron[2] = &Necker[5+8]; Necker[7+8].inp_wgt[2] = Excitation;
	Necker[7+8].inp_neuron[3] = &Necker[6+8]; Necker[7+8].inp_wgt[3] = Excitation;

}

int main(int argc, char **argv)
{
	/* SEED and loadingtime are required in order to get random numbers */
	unsigned int SEED;
	long loadingtime;

	/* initialize of random number sequence */
	if (argc == 2) {
		if ((SEED = atoi(argv[1])) == 0) {
			fprintf(stderr, "Invalid arg :%s\n",argv[1]);
			exit (EXIT_FAILURE);
		}
	} else {
		SEED = (unsigned)time(&loadingtime); 
	}
	srand(SEED);
	printf("### Randomd SEED: %d\n", SEED);

	/* Make neural network for Necker Cube */
	Make_Necker_Cube_Network();
	main_loop();
	return EXIT_SUCCESS;
}