/* coupled neurons */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>  /* time.h is required for random number generator */

#define FROM   0.0
#define TO   100.0
#define STEP   0.01

#define sigmoid_slope -1.0

double INPUT=0.1;
double INITIAL_VALUE=10.0;
double INIT_SELFWGT=1.0;

/* SEED and loadingtime are required in order to get random numbers */
unsigned int SEED;
long loadingtime;

/* Definition of a neuron
 * IMPORTANT!!!!  Among of the member of the following structure,
 *  Ninp strands for a number of input from neurons.
 */
struct NEURON {
    double output, old_output;
    double thres;
    int    Ninp;
    double *inp_wgt;
    struct NEURON **inp_neuron;
};
typedef struct NEURON neuron;

/* Each neuron can be defined as global variables */
neuron a, b, to_a_neuron, to_b_neuron;

/* 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 = drand2();
    a->output = drand2();
    a->Ninp = N;
    a->inp_wgt = (double *)malloc(sizeof(double) * (N+1));
    for (i=0; i<=N; i++) {
	a->inp_wgt[i] = drand2();
    }
    a->inp_neuron = (neuron **)malloc(sizeof(neuron *) * (N+1));
    /* self feedback */
		/* a->inp_neuron[0] = a; */
}

double output_f(double value)
{
    double ret;
    ret = 1.0 / ( 1.0 + exp(sigmoid_slope * value) );
    /*
    return ret;
    */
    return value;
}

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

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

void update_neuron(neuron *a)
{
    a->old_output = a->output;
}

int main_loop(void)
{
    double wrk;
    double t = FROM;
    while ( t <= TO ){
	calc_neuron(&a); 	calc_neuron(&b); 
	update_neuron(&a); 	update_neuron(&b); 
	printf("%5.2f %7.4f %7.4f\n", t, a.output, b.output);
	t += STEP;
    }
    return 0;
}

int main(int argc, char **argv)
{
    if ( argc != 3 ) {
	fprintf(stderr,"Usage: %s init_a init_b\n",argv[0]);
	exit( 0 );
    }

    /* initialize of random number sequence */
    SEED = (unsigned)time(&loadingtime); srand(SEED);

    initialize_neuron(&a, 2);
    initialize_neuron(&b, 2);

    initialize_neuron(&to_a_neuron, 0);
    initialize_neuron(&to_b_neuron, 0);

    a.inp_neuron[0] = &to_a_neuron;
    b.inp_neuron[0] = &to_b_neuron;
    
    a.inp_neuron[1] = &b;
    b.inp_neuron[1] = &a;

    /*
    to_a_neuron.output = to_a_neuron.old_output = -3.0;
    to_b_neuron.output = to_b_neuron.old_output = 4.0;
    */
    to_a_neuron.output = to_a_neuron.old_output = atof(argv[1]);
    to_b_neuron.output = to_b_neuron.old_output = atof(argv[2]);

    /*
    a.output = 1.0; 
    b.output = -1.0;
    */

    a.inp_wgt[0] = b.inp_wgt[0] = 1.0;
    a.inp_wgt[1] = 0.0;
    b.inp_wgt[1] = 0.0;
    a.thres = b.thres = 0.0;

#ifdef DEBUG
    fprintf(stderr,"a.output=%7.4f, a.inp_neuron[0]->output=%g, to_a_neuron.output=%g, a.inp_neuron[1]->output=%g\n", 
	   a.output, a.inp_neuron[0]->output, to_a_neuron.output,
	   a.inp_neuron[1]->output);
    fprintf(stderr,"b.output=%7.4f, b.inp_neuron[0]->output=%g, to_b_neuron.output=%g, b.inp_neuron[1]->output=%g\n", 
	   b.output, b.inp_neuron[0]->output, to_b_neuron.output,
	   b.inp_neuron[1]->output);

    fprintf(stderr,"a:thres=%g, inp_wgt[0]=%g, inp_wgt[1]=%g, init=%g\n",
	   a.thres, a.inp_wgt[0], a.inp_wgt[1], a.output);
    fprintf(stderr,"b:thres=%g, inp_wgt[0]=%g, inp_wgt[1]=%g, init=%g\n",
	   b.thres, b.inp_wgt[0], b.inp_wgt[1], b.output);
#endif 
    
    return main_loop();
}