/* boltzoutput.c: Generate outputs for the already trained boltzmann machine
*/
#include <sys/time.h>
#include <stdio.h>
#include <math.h>

/* Maximum random number used for the purposes generating a 
   double-precision floating point random number */
#define MAXRAND (1<<16)
#define NOISE_OFF2ON 0.05
#define NOISE_ON2OFF 0.05

/* Storage for the machine. All data relevant for the machine are kept in 
   this structure */
struct machine
{	/* Machine data */
	double ***weights; /* weights[layer i][node in layer i]
					      [node in layer i+1] */
	double ***internal; /* intra-layer weights. 
				internal[layer i][from node1 in layer i]
						[to node2 in layer i] 
				The matrix is symmetrical across it's second
				and third dimensions! Therefore, only
				those weights where node1>node2 are stored. 
				The symmetric weight is equal to that */
	int *internal_exist;  /* Do intralayer connections exist for each
				 layer: internal_exist[layer] */
	int *num_nodes;	/* Number of nodes in each layer: num_nodes[layer] */
	int layers;	/* Number of layers in the machine. Should be 3 */
	int **states;	/* State of the machine: states[layer][node] */
	int *clamped;	/* Is the i-th layer clamped, i.e. fixed */

	/* Input storage */
	int num_inputs; /* Number of inputs */
	int **inputs;	/* The inputs the machine will be trained on: 
				inputs[input#][dimension] */
	int **outputs;	/* The outputs that correspond to the training
				inputs. outputs[output#][dimension]

	/* Training algorithm data */
	int max_iter;		/* Maximum iterations */
	double eta;		/* learning rate */
	double converge;	/* Convergence criterion. */
	double t0, beta, tmin;	/* annealing start and end temperatures */
	
	/* Mean field annealing data */
	int mean_field;		/* Is the annealing mean_field? */
	double **mean_states; 	/* mean values of states */
	double **pon, **poff; 	/* prob. of states being on or off */

	/* Training algorithm probabilities */
	/* probabilities gathered in phases 1 and 2 (p1 and p2) on 
	   the weights and the internal weights 
	*/
	double ***p1weights, ***p1internal, 
	       ***p2weights, ***p2internal;	
};


/* common functions */
int get_1bit(), allocate_weights();
double get_rand(), energy_change();
struct machine *get_machine();
void init_random(), deallocate_machine();
void flip_state(), noisy_flip_state(), select_node(), init_state();
int get_num_changes(), accept_change();
void anneal(), anneal_1_step();
void mean_field_anneal();
double get_mean_input();
void init_exp_lookup();
double exp_lookup();


/* Common functions for training and getting output from a boltzmann machine
*/


#define ON (1)
#define OFF (-1)

/* Initialize random generator with number of microseconds since
   January 1, 1970, which is almost a random value in itself */
void initrandom()
{
	struct timeval tp;
	struct timezone tzp;
	if( gettimeofday(&tp,&tzp) == -1 )
	{       printf("Could not initialize random number generator\n");
		_exit(0);
	}
	srandom((int)tp.tv_usec);
}


/* Returns a double precision random number between low and high  */
double get_rand(low, high)
double low, high;
{	double x = random() % MAXRAND / (double)MAXRAND;
	return( x*(high-low) + low );
}


/* Returns a 1-bit number */
int get_1bit()
{	return( random() % 2 );
}


/* deallocates a weights matrix */
deallocate_weight_matrix(m,p)
double ***m;
struct machine *p;
{	int i,j;

	if( m )
	{   for( i=0; i<p->layers-1; i++ )
	    {	if( m[i] )
		{   for( j=0; j<p->num_nodes[i]; j++)
			if( m[i][j] )
			    free( m[i][j] );
		    free(m[i]);
		}
	    }
	    free(m);
	}
}

/* deallocates an internal weights matrix */
deallocate_internal_weight_matrix(m,p)
double ***m;
struct machine *p;
{	int i,j;

	if( m )
	{   for( i=0; i<p->layers; i++ )
	    {	if( m[i] )
		{   for( j=0; j<p->num_nodes[i]; j++)
			if( m[i][j] )
			    free( m[i][j] );
		    free(m[i]);
		}
	    }
	    free(m);
	}
}

/* deallocates one of the structures associated with mean field annealing */
deallocate_mean_field( m, p)
double **m;
struct machine *p;
{	int i;

	if( m )
	{   for( i=0; i<p->layers; i++ )
		if( m[i] )
		    free( m[i] );
	    free(m);
	}
}

/* deallocates the structure and all the allocated arrays within */
void deallocate_machine(p)
struct machine *p;
{	int i,j;

	if( p->num_nodes )
		free(p->num_nodes);

	deallocate_weight_matrix(p->weights, p);

	if( p->internal_exist )
		free( p->internal_exist );

	deallocate_internal_weight_matrix(p->internal, p);

	if( p->states )
	{   for( i=0; i<p->layers; i++ )
		if( p->states[i] )
		    free( p->states[i] );
	    free(p->states);
	}

	if( p->clamped )
	    free( p->clamped );

	if( p->inputs )
	{   for( i=0; i<p->num_inputs; i++ )
		if( p->inputs[i] )
		    free(p->inputs[i]);
	    free(p->inputs);
	}

	if( p->outputs )
	{   for( i=0; i<p->num_inputs; i++ )
		if( p->outputs[i] )
		    free(p->outputs[i]);
	    free(p->outputs);
	}

	deallocate_weight_matrix(p->p1weights, p);
	deallocate_weight_matrix(p->p2weights, p);
	deallocate_internal_weight_matrix(p->p1internal, p);
	deallocate_internal_weight_matrix(p->p2internal, p);

	deallocate_mean_field( p->mean_states, p );
	deallocate_mean_field( p->pon, p );
	deallocate_mean_field( p->poff, p);

	free(p);
}


/* Allocates a structure for the machine, and initializes all dynamic
   strutures to NULL */
struct machine *get_machine()
{	struct machine *p;

	p = (struct machine *)malloc(sizeof(struct machine));
	p->weights = p->internal = NULL;
	p->internal_exist = NULL;
	p->num_nodes = NULL;
	p->clamped = NULL;
	p->inputs = p->outputs = NULL;
	p->p1weights = p->p2weights = p->p1internal = p->p2internal = NULL;
	p->mean_states = p->pon = p->poff = NULL;
	return(p);
}


/* Allocate weight matrix (according to data in *p), and return 1 if 
   successfull. If not, return 0. This function will be used for allocating
   the weight matrix, as well as the 2 probability matrices associated
   with the weight matrix
*/
allocate_weight_matrix(pw,p)
double ****pw;
struct machine *p;
{	double ***w;
	int i,j;

	/* Allocate weights matrix */
	*pw = (double ***)malloc((p->layers-1)*sizeof(double **));
	if( *pw == NULL )
		return(0);
	w = *pw;
	for( i=0; i<p->layers-1; i++ )
		w[i] = NULL;	/* So that deallocation of partially 
					   allocated weights is possible */
	for( i=0; i<p->layers-1; i++ )
	{    w[i] = (double **)malloc(p->num_nodes[i]*sizeof(double *));
	     if( w[i] == NULL )
		return(0);
	     for( j=0; j<p->num_nodes[i]; j++ )
		w[i][j] = NULL; 
	     for( j=0; j<p->num_nodes[i]; j++ )
	     {  w[i][j] = (double *)malloc(p->num_nodes[i+1]*sizeof(double)); 
		if( w[i][j] == NULL )
			return(0);
	     }
	}

	/* Success */
	return(1);
}

/* Allocate internal weight matrix (according to data in *p), and return 
   1 if successfull. If not, return 0. This function will be used for 
   allocating the internal matrix, as well as the 2 probability matrices 
   associated with the internal matrix
*/
allocate_internal_weight_matrix(pw,p)
double ****pw;
struct machine *p;
{	double ***w;
	int i,j;

	/* Allocate internal weights matrix */
	*pw = (double ***)malloc(p->layers*sizeof(double **));
	if( *pw == NULL )
		return(0);
	w = *pw;
	for( i=0; i<p->layers; i++ )
		w[i] = NULL;	/* So that deallocation of partially 
					   allocated weights is possible */
	for( i=0; i<p->layers; i++ )
	  if( p->internal_exist[i] )		/* allocated only if needed */
	  {  w[i] = (double **)malloc(p->num_nodes[i]*sizeof(double *));
	     if( w[i] == NULL )
		return(0);
	     for( j=0; j<p->num_nodes[i]; j++ )
		w[i][j] = NULL; 
	     for( j=0; j<p->num_nodes[i]; j++ )
	     {  w[i][j] = (double *)malloc(p->num_nodes[i]*sizeof(double)); 
		if( w[i][j] == NULL )
			return(0);
	     }
	  }

	/* Success */
	return(1);
}


/* Allocates the 3 arrays associated with the mean field annealing */
allocate_mean_field(p)
struct machine *p;
{	int i;

	/* Allocate the mean states matrix */
	p->mean_states = (double **)malloc(p->layers*sizeof(double *));
	if( p->mean_states == NULL )
		return(0);
	for( i=0; i<p->layers; i++ )
		p->mean_states[i] = NULL;
	for( i=0; i<p->layers; i++ )
	{    p->mean_states[i] = (double *)malloc(p->num_nodes[i]*sizeof(double));
	     if( p->mean_states[i] == NULL )
		return(0);
	}

	/* Allocate the pon matrix */
	p->pon = (double **)malloc(p->layers*sizeof(double *));
	if( p->pon == NULL )
		return(0);
	for( i=0; i<p->layers; i++ )
		p->pon[i] = NULL;
	for( i=0; i<p->layers; i++ )
	{    p->pon[i] = (double *)malloc(p->num_nodes[i]*sizeof(double));
	     if( p->pon[i] == NULL )
		return(0);
	}

	/* Allocate the poff matrix */
	p->poff = (double **)malloc(p->layers*sizeof(double *));
	if( p->poff == NULL )
		return(0);
	for( i=0; i<p->layers; i++ )
		p->poff[i] = NULL;
	for( i=0; i<p->layers; i++ )
	{    p->poff[i] = (double *)malloc(p->num_nodes[i]*sizeof(double));
	     if( p->poff[i] == NULL )
		return(0);
	}

	/* Success */
	return(1);
}


/* Initializes the arrays in the structure 'machine'. All it's parameters
   are already in the structure, such as: #layers, #nodes in each layer 
   and the internal_exist array. If any allocation does not succeed, 
   returns 0 (false)
*/
int allocate_weights(p)
struct machine *p;
{
	int i,j,k;

	/* Allocate the weight matrix */
	if( allocate_weight_matrix(&(p->weights), p) == 0)
		return(0);

	/* Allocate the internal weights matrix */
	if( allocate_internal_weight_matrix(&(p->internal), p) == 0)
		return(0);

	/* Allocate the states matrix */
	p->states = (int **)malloc(p->layers*sizeof(int *));
	if( p->states == NULL )
		return(0);
	for( i=0; i<p->layers; i++ )
		p->states[i] = NULL;
	for( i=0; i<p->layers; i++ )
	{    p->states[i] = (int *)malloc(p->num_nodes[i]*sizeof(int));
	     if( p->states[i] == NULL )
		return(0);
	}

	/* Allocate the 'clamped' array */
	p->clamped = (int *)malloc(p->layers * sizeof(int));
	if( p->clamped == NULL )
		return(0);

	/* Allocate the inputs matrix */
	p->inputs = (int **)malloc(p->num_inputs*sizeof(int *));
	if( p->inputs == NULL )
		return(0);
	for( i=0; i<p->num_inputs; i++ )
		p->inputs[i] = NULL;
	for( i=0; i<p->num_inputs; i++ )
	{    p->inputs[i] = (int *)malloc(p->num_nodes[0]*sizeof(int));
	     if( p->inputs[i] == NULL )
		return(0);
	}

	/* Allocate the outputs matrix. There are as many outputs 
	   as there are inputs in the machine */
	p->outputs = (int **)malloc(p->num_inputs*sizeof(int *));
	if( p->outputs == NULL )
		return(0);
	for( i=0; i<p->num_inputs; i++ )
		p->outputs[i] = NULL;
	for( i=0; i<p->num_inputs; i++ )
	{    p->outputs[i] = (int *)malloc(p->num_nodes[p->layers-1]*sizeof(int));
	     if( p->outputs[i] == NULL )
		return(0);
	}

	/* Allocate probability matrices */
	if( allocate_weight_matrix(&(p->p1weights),p) == 0 )
		return(0);
	if( allocate_weight_matrix(&(p->p2weights),p) == 0 )
		return(0);
	if( allocate_internal_weight_matrix(&(p->p1internal),p) == 0 )
		return(0);
	if( allocate_internal_weight_matrix(&(p->p2internal),p) == 0 )
		return(0);

	/* Success! */
	return(1);
}


/* flips the state of node j in layer i */
void flip_state(p,i,j)
struct machine *p;
int i,j;
{
	p->states[i][j] = (p->states[i][j] == ON) ? OFF : ON;
}


/* Flips the state of p->states[i][j] with a certain probability. 
   probabilities defined in the header file */
void noisy_flip_state(p,i,j)
struct machine *p;
int i,j;
{
	if( p->states[i][j] == ON && get_rand(0.0, 1.0) < NOISE_ON2OFF )
		p->states[i][j] = OFF;
	else if( p->states[i][j] == OFF && get_rand(0.0, 1.0) < NOISE_OFF2ON )
		p->states[i][j] = ON;
}


/* randomly selects one node. Selections should be based on clamped 
  inputs, but that is the second stage. 
*/
void select_node(p, layer, node)
struct machine *p;
int *layer, *node;
{
	int sum=0, i;
	int num;
	int selected = 0;

	/* Randomly select a number between 0 and the total number of nodes, 
	   not counting the clamped nodes */
	for( i=0; i<p->layers; i++)
		sum += (1-p->clamped[i])*p->num_nodes[i];

	/* Keep trying until a node that is not in one of the clamped layers
	   is selected */
	num = (int) get_rand(0.0, (double)(sum) - 0.001);
	
	/* find the corresponding layer and node. Nodes are 'numbered' 
	   starting from 0 (first node in first layer), to sum-1 (last 
	   node in last layer. Only unclamped nodes are 'valid' */
	for( i=0; i<p->layers && num >= 0 ; i++ )
	    if( !p->clamped[i] )
	    {	if( num < p->num_nodes[i] )
		{	*layer = i;
			*node = num;
		}
		num -= p->num_nodes[i];
	    }
}


/* Compute the change of energy if states[layer][node] were flipped */
double energy_change(p, layer, node)
struct machine *p;
int layer, node;
{	
	int old,new;
	int i,j;
	double sum = 0;
	/* When state is OFF or ON, the change of energy is
	   (new_state - old_state) * sum(w*x, for all weights leading 
			into the node, except for the self-excitation weight) 
	   for an energy defined as: -sum[w(i,j)*x(i)*x(j)]
	*/
	/* Get the old state, and flip it */
	old = p->states[layer][node];
	flip_state(p, layer, node);
	new = p->states[layer][node];

	/* 3 types of weights lead into node (layer,node): 
	   - weights from (layer-1) to (layer), if such exist
	   - weights from (layer) to (layer+1), if such exist
	   - weights from (layer) to (layer), if such exist
	*/
	/* Connections from prev. layer */
	if( layer-1 >= 0 )
	    for( j=0; j<p->num_nodes[layer-1]; j++ )
		sum += p->weights[layer-1][j][node] * p->states[layer-1][j];

	/* Connections to next layer */
	if( layer < p->layers-1 )
	    for( j=0; j<p->num_nodes[layer+1]; j++)
		sum += p->weights[layer][node][j] * p->states[layer+1][j];

	/* Connections in the same layer. The internal connections are 
	   stored in internal[layer][a][b] iff a>b. If b>a, then they are in 
	   internal[layer][b][a]. The reason is that internal is a symmetric 
	   matrix 
	*/
	if( p->internal_exist[layer] )
	{   for( j=0; j<p->num_nodes[layer]; j++ )
		if( j < node )
		    sum += p->internal[layer][node][j] * p->states[layer][j];
		else if( j > node )
		    sum += p->internal[layer][j][node] * p->states[layer][j];
	}
	/* Restore old state */
	p->states[layer][node] = old;
/*
printf("Sum=%lf\n",sum);

printf("New %d Old %d\n",new,old);
printf("New-old %d\n",new-old);
printf("Sum of weights and states for [%d][%d] is %lf\n",layer,node,sum);
printf("Returning %lf\n",(new-old)*sum);
*/
	return((new-old)*sum);
}


/* Initialize the state of the machine by setting all nodes to random 
   states */
void initstate(p)
struct machine *p;
{	int i,j;

	/* All nodes in state ON */
	for( i=1; i<p->layers; i++)
	    if( !p->clamped[i] )
		for( j=0; j<p->num_nodes[i]; j++)
		    p->states[i][j] = ON;
	/* Randomly flip state of  all nodes in layer diff. than input */
	for( i=1; i<p->layers; i++)
	    if( !p->clamped[i] )
		for( j=0; j<p->num_nodes[i]; j++)
		    if( get_1bit() )
			flip_state(p,i,j);
}


/* Which number of energy changes does the machine perform for a given 
   temperature */
int get_num_changes(t,p)
double t;
struct machine *p;
{	double num = p->t0 / t;
	int i;
	double total=0;

	for( i=0; i<p->layers; i++)
		total += (1-p->clamped[i]) * num * p->num_nodes[i];

	return((int)total);
}

/* Is a change of dE acceptable at temperature t? */
accept_change(dE,t)
double dE,t;
{	double prob, rand;

	/* Always accept changes that increase the energy */
	if( dE > 0 )
		return( 1 );

	/* If the change decreases the energy, accept it with a
	   certain probability */

	prob = 1 / (1 + exp_lookup(-dE/t) );
	rand = get_rand(0.0, 1.0);
	return( rand < prob );

}


/* Simulated annealing over machine p. The temperature is constantly 
   decreasing over a set of values. For each temperature a set of state 
   changes are performed on randomly selected (non clamped!) nodes. Each 
   state change is selected with a certain probability. If it decreases 
   the total energy, it is selected, if not it is selected with a 
   probability that is a function of the temperature
*/
void anneal(p)
struct machine *p;
{	double temp = p->t0;
	int i, n;
	int node, layer;
	double dE;

	/* Vary the temperature */ 
	while( temp >= p->tmin )
	{   /* For each temperature, perform a number of operations which 
	       is a function of the temperature of annealing. */
	    n = get_num_changes(temp, p);
	    for( i=0; i<n; i++ )
	    {	select_node(p, &layer, &node);
		dE = energy_change(p, layer, node);
		if( accept_change(dE, temp) )
		    flip_state(p, layer, node);
	    }
	    temp *= p->beta;
	}
}


/* Change the value for only 1 node, on temperature t. On this temperature 
   accept the change if it increases the energy, and accept it with
   some probability, if it decreases it. Probability depends on t */
void anneal_1_step(p,t)
struct machine *p;
double t;
{	int node, layer;
	double dE;

	select_node(p, &layer, &node);
	dE = energy_change(p, layer, node);
	if( accept_change(dE, t) )
		flip_state(p, layer, node);
}


/* Lookup table for values of exp(x). x starts from 0 and increases in 
   steps of 0.02 */
double evalues[750];

/* Initialize lookup table */
void init_exp_lookup()
{	double n;
	int i;

	for( i=0; i<750; i++ )
	{	n = (double) i / 50;
		evalues[i] = exp(n);
	}
}

/* Lookup the value for exp(n) */
double exp_lookup(n)
double n;
{	int i;

	i = n * 50;
	return( (i>499) ? evalues[499] : evalues[i] );
}


/* Get the net input to node [layer][node]. Net input is calculated 
   as a function of the mean states of the surrounding nodes
*/
double get_mean_input(p, layer, node)
struct machine *p;
int layer, node;
{	double sum = 0;
	int j;      

	/* Connections from prev. layer */
	if( layer-1 >= 0 )
	    for( j=0; j<p->num_nodes[layer-1]; j++ )
		sum += p->weights[layer-1][j][node] * 
					p->mean_states[layer-1][j];

	/* Connections to next layer */
	if( layer < p->layers-1 )
	    for( j=0; j<p->num_nodes[layer+1]; j++)
		sum += p->weights[layer][node][j] * 
					p->mean_states[layer+1][j];
	/* Connections in the same layer */
	if( p->internal_exist[layer] )
	{   for( j=0; j<p->num_nodes[layer]; j++ )
		if( j < node )
		    sum += p->internal[layer][node][j] * 
					p->mean_states[layer][j];
		else if( j > node )
		    sum += p->internal[layer][j][node] * 
					p->mean_states[layer][j];
	}

	return(sum);
}


/* Gets mean state of a node, depending on the sum and temperature. 
   Is a mapping of the tanh(sum/temp) from the (-1,+1) to the (OFF,ON) 
   interval. 
 */
double get_mean_state( sum, temp)
  double sum, temp;
{	double x;

	x = tanh(sum/temp);		/* Compute tanh */
/*	x = (x+1);			/* Map to (0,2) 
	x = x * (ON-OFF) / 2.0;		/* Map to (0,(ON-OFF)) 
	x = x + OFF;			/* Map to (OFF,ON) 
*/
	return(x);
}


/* mean field annealing algorithm. Perform mean field annealing for each 
   temperature. for now: just one pass to determine the average values of 
   the nodes 
*/
void mean_field_anneal(p)   
struct machine *p;
{	double temp = p->t0;
	double sum;
	int layer, node, i, j;
	int repetitions=2;

	/* assign real values to average values */
	for( i=0; i<p->layers; i++)
	    for( j=0; j<p->num_nodes[i]; j++)
		p->mean_states[i][j] = p->states[i][j];

	/* For all temperatures */
	while( temp >= p->tmin )
	{
	  /* repeat a few times */
	  for( i=0; i<repetitions; i++)
	    /* For each node */
	    for( layer=0; layer<p->layers; layer++)
		if( !p->clamped[layer] )
		    for( node=0; node<p->num_nodes[layer]; node++)
		    {	sum = get_mean_input(p, layer, node);
			p->mean_states[layer][node] = get_mean_state(sum, temp); 
		    }

	    temp *= p->beta;
	}
	/* At the end of the "annealing", find the probabilities, 
	   and update the node states accordingly!  */
	for( layer=0; layer<p->layers; layer++)
	    for( node=0; node<p->num_nodes[layer]; node++)
	    {   p->pon[layer][node] = (p->mean_states[layer][node] - OFF);
		p->pon[layer][node] /= (double) (ON-OFF);

		p->poff[layer][node] = 1 - p->pon[layer][node];

		if( p->pon[layer][node] > p->poff[layer][node] )
		    p->states[layer][node] = ON; 
		else
		    p->states[layer][node] = OFF;
	    }
}



void get_user_parameters();
void boltzoutput(), openfiles();
struct machine *readinput();
void writeoutput();
void hopfield_style(), anneal_style();


/* argv[1] - input weights file
   argv[2] - input data file
*/
void main(argc, argv)
int argc;
char *argv[];
{	int typ;
	int i=0;
	char weightfile[50], testfile[50];
	char type;

	if( argc > 1 && argv[1][0] == '-' )
	{	typ = argv[1][1];
		i=1;
	}		
	else
		typ = 'a';

	if( argc - i == 1 )		/* Query user for values */
	{	get_user_parameters(weightfile, testfile, &type);
		typ = type;
		boltzoutput(typ, weightfile, testfile);
	}
	else
		boltzoutput(typ, argv[i+1], argv[i+2]);
}


/* Prompts user for names of files with weights and test data, as well as the
   desired output-generation algorithm 
*/
void get_user_parameters( wfile, tfile, type)
char *wfile, *tfile, *type;
{
	printf("Enter the name of the weight file: ");
	scanf("%s",wfile);
	printf("Enter the name of the test data file: ");
	scanf("%s",tfile);
	printf("Select output generation algorithm: \n");
	printf("Type \n h for hopfield \n a for annealing \n");
	printf(" m for mean-field annealing\n");
	scanf("\n%c",type);
}


/* Opens the files. If an error occurs in opening the files, the program
   exits */
void openfiles(infile, datafile, in, indata)
char *infile, *datafile;
FILE **in, **indata;
{
	if( (*in = fopen(infile,"rt")) == NULL )
	{	printf("Cannot open file %s for reading\n",infile);
		_exit(0);
	}
	if( (*indata = fopen(datafile,"rt")) == NULL )
	{	printf("Cannot open file %s for reading\n",datafile);
		_exit(0);
	}
}


/* Presents each input to a boltzmann machine, and gets an output 
   using the hopfield algorithm.
*/
hopfield_output(p)
struct machine *p;
{
	int i;

	initrandom();
	init_exp_lookup();
	for( i=0; i<p->num_inputs; i++ )
	{   hopfield_style(p, p->inputs[i] );
	    writeoutput(p);
	}
}


/* Generates outputs from the boltzmann machine by a hopfield-style algorithm. 
   At start, the input nodes receive the input pattern, and all other nodes 
   receive random values. State of the machine is changed node by node (by minimizing
   the energy function), until no further changes can be made. 
*/
void hopfield_style(p,inputs)
struct machine *p;
int *inputs;
{
	int i,j;
	int layer,node; 
	int inlayer = 0;
	int nochanges = 0;

	/* Initialize to random state */
	initstate(p);

	/* Copy inputs to the nodes in the input layer (0-th layer) */
	for( i=0; i<p->num_nodes[0]; i++ )
	    p->states[0][i] = inputs[i];
	/* Clamp only input nodes */
	p->clamped[inlayer] = 1;

	/* While changes are possible, select 1 node, and see whether it can
	   be changed by examining the difference in the energy f-n
	*/
	i=0;
	while( nochanges < 20 )
	{   select_node(p, &layer, &node);
	    if( energy_change(p, layer, node) > 0 )
	    {	flip_state(p, layer, node);
		nochanges = 0;
	    }
	    else nochanges++;
	    i++;
	}
	/* Release input nodes */
	p->clamped[inlayer] = 0;
}


/* Presents each output to a boltzmann machine and gets an output by setting a
   full annealing process
*/
anneal_output(p)
struct machine *p;
{	int i;

	initrandom();
	init_exp_lookup();
	for( i=0; i<p->num_inputs; i++ )
	{   anneal_style(p, p->inputs[i] );
	    writeoutput(p);
	}
}


/* Generates output from the boltzmann machine using an annealing algorithm. 
   At start an input pattern is presented to the machine. The states of the machine
   are changed according to the annealing algorithm, i.e. changes that increase 
   the energy are allowed, and reverse changes are allowed only with certain 
   probability. 
*/
void anneal_style(p,inputs)
struct machine *p;
int *inputs;
{	int i;
	int inlayer = 0;

	/* Initialize to random state */
	initstate(p);

	/* Copy inputs to the nodes in the input layer (0-th layer) */
	for( i=0; i<p->num_nodes[0]; i++ )
	    p->states[0][i] = inputs[i];
	/* Clamp only input nodes */
	p->clamped[inlayer] = 1;

	if( p->mean_field )
	{	p->t0 = 40.0;
		p->tmin = 0.5;
		p->beta = 0.9;
		mean_field_anneal(p);
	}
	else
	{
		p->t0 = 20.0;
		p->tmin = 0.1;
		p->beta = 0.95;
		anneal(p);
	}

	/* Release input nodes */
	p->clamped[inlayer] = 0;
}


/* Reads the machine state (the structure, the weights) from the 'in' file. On 
   unsuccessfull read, a NULL pointer is returned. On successfull read, a pointer
   towards a structure 'machine' containing all the info is returned 
*/
struct machine *readinput(in, indata)
FILE *in, *indata;
{	struct machine *p;
	int i,j,k;

	/* Initialize the machine */
	if( (p = get_machine()) == NULL )
		return(NULL);
	        
	/* Read file and allocate arrays. On unsuccesfull allocation,
	   deallocate all and return a NULL pointer, signifying non-success
	*/
	fscanf(in, "%d\n", &(p->layers));
	if( (p->num_nodes = (int *)malloc(p->layers * sizeof(int))) == NULL)
	{	deallocate_machine(p);
		return(NULL);
	}

	/* Read number of nodes in each layer */
	for( i=0; i<p->layers; i++ )
		fscanf(in, "%d ",p->num_nodes+i);
	fscanf(in,"\n");


	/* Read the existence matrix */
	p->internal_exist = (int *)malloc(p->layers*sizeof(int));
	if( p->internal_exist == NULL )
	{	deallocate_machine(p);
		return(NULL);
	}
	for( i=0; i<p->layers; i++ )
		fscanf(in, "%d ", p->internal_exist+i);
	fscanf(in,"\n");

	/* Read the number of inputs */
	fscanf(indata, "%d\n", &(p->num_inputs));

	/* allocate space for weights, internal weights, the state of the 
	   machine, the inputs etc. */
	if( !allocate_weights(p) )
	{	deallocate_machine(p);   
		return(NULL);
	}

	/* Just in case: allocate mean_field matrices */
	if( !allocate_mean_field(p))
	{	deallocate_machine(p);
		return(NULL);
	}

	/* Read in the weights */
	for( i=0; i<p->layers-1; i++)
	    for( j=0; j<p->num_nodes[i]; j++)
		for( k=0; k<p->num_nodes[i+1]; k++)
		    fscanf(in, "%lf\n", &(p->weights[i][j][k]));
	/* Read in the internal weights */
	for( i=0; i<p->layers; i++)
	    if( p->internal_exist[i] )
		for( j=0; j<p->num_nodes[i]; j++)
		    for( k=0; k<j; k++)
			fscanf(in, "%lf\n", &(p->internal[i][j][k]));

	/* Read in the inputs */
	for( i=0; i<p->num_inputs; i++ )
	{   for( j=0; j<p->num_nodes[0]; j++ )
		fscanf(indata, "%d ", p->inputs[i] + j );
	    fscanf(indata, "\n");
	}

	/* Mark all layers as unclamped */
	for( i=0; i<p->layers; i++ )
	    p->clamped[i] = 0;

	/* That's all the reaing */
	return(p);
}


/* should write out the state of the machine after each output is generated */
void writeoutput(p)
struct machine *p;
{
	int i,j,k;

printf("Node states: \n");
for( i=0; i<p->layers; i++)
{	printf("Layer %d: ",i);
	for( j=0; j<p->num_nodes[i]; j++)
		printf("%2d ",p->states[i][j]);
	printf("\n");
}

}


/* Generate outputs for the boltzmann machine whose weights are stored in file
   'infile'
*/
void boltzoutput(typ, infile, datafile)
int typ;
char *infile, *datafile;
{	FILE *in, *indata;
	struct machine *p;

	openfiles(infile, datafile, &in, &indata);
	if( (p=readinput(in, indata)) == NULL)
	{	printf("Not enough memory\n");
		_exit(0);
	}

	switch(typ) 
	{	case 'a': 
			p->mean_field = 0;
			anneal_output(p);
			break;
		case 'm': 
			p->mean_field = 1;
			anneal_output(p);
			break;
		case 'h': 
			hopfield_output(p);
			break;
		default: 
			printf("Cannot generate output with the -%c option\n",typ);
			break;
	}
}