// $Id$
//
//    File: nnann.cc
// Created: Fri Mar 18 19:47:27 EDT 2011
// Creator: davidl (on Darwin Amelia.local 9.8.0 i386)
//

#include <nnann.h>
#include <nnann_function_sigmoid.h>

using namespace std;


//---------------------------------
// nnann    (Constructor)
//---------------------------------
nnann::nnann():activation_function(NULL),initialized(false)
{
	// A single nnann_neuron is used for the bias. The neuron always emits an
	// ouput of 1.0 allowing the weights to represent the real bias in every 
	// neurons input. The bias neuron is included in the connection layers, but
	// not in the neuron layers themselves. (See RelinkLayers.)
	bias.SetOutput(1.0);
}

//---------------------------------
// nnann    (Constructor)
//---------------------------------
nnann::nnann(std::string filename):activation_function(NULL),initialized(false)
{
	ReadFile(filename);
}

//---------------------------------
// ~nnann    (Destructor)
//---------------------------------
nnann::~nnann()
{

}

//---------------------------------
// AddLayer
//---------------------------------
void nnann::AddLayer(unsigned int Nneurons, nnann_function *activation_function)
{
	// Create a layer with specified number of neurons and the specified activation function
	nnann_layer layer(Nneurons);
	for(unsigned int i=0; i<layer.size(); i++)layer[i].SetActivationFunction(activation_function);

	// Add this layer to the ANN
	layers.push_back(layer);

	// Relink all of the layers since the above push_back could cause memory reallocation
	if(layers.size()>1)RelinkLayers();	
}

//---------------------------------
// BackPropagate
//---------------------------------
double nnann::BackPropagate(const std::vector<double>& target_vals)
{
	/// Use the specified target_vals to calculate the error on the output layer and
	/// propagate it back through the ANN, updating the weights. The mean squared
	/// value of the output errors is returned (i.e. the "MS" not the "RMS").

	// Dimension check
	if(layers.size()<2)throw nnann_exception("Can't back propagate. Network has less than 2 layers!");
	nnann_layer &output_layer = layers[layers.size()-1];
	if(output_layer.size()!=target_vals.size())throw nnann_exception("Can't back propagate. target vector size does not match network's output layer size!");

	// Learning rate (see comment below)
	double eta = 0.03;

	// Clear the output errors for all neurons in the ANN
	for(unsigned int ilayer=0; ilayer<layers.size(); ilayer++)
		for(unsigned int ineuron=0; ineuron<layers[ilayer].size(); ineuron++)layers[ilayer][ineuron].ResetOutputError();

	// Set output layer errors to the difference between them and the target values
	// Keep track of sum of squares of errs so we can return the mean squared value.
	double sum_squared_error = 0.0;
	for(unsigned int ineuron=0; ineuron<output_layer.size(); ineuron++){
		double err = target_vals[ineuron] - output_layer[ineuron].GetOutput();
		output_layer[ineuron].SetOutputError(err);
		sum_squared_error += err*err;
	}

	// Loop over all connection layers, propagating the errors backward.
	// n.b. We don't calculate errors on the input layer of neruons as these
	// are given and so can't be adjusted by the network.
	for(unsigned int ilayer=(weights_layers.size()-1); ilayer>0; ilayer--){
	
		// Loop over connections in this connection layer, adding to the output error of the previous layer
		nnann_weight_layer &weights = weights_layers[ilayer];
		for(unsigned int iweight=0; iweight<weights.size(); iweight++){
			weights[iweight].BackPropagate();
		}
	
		// Loop over neurons in the layer we just propagated errors to and apply 
		// the activation function making their input errors valid
		vector<nnann_neuron> &layer = layers[ilayer];
		for(unsigned int ineuron=0; ineuron<layer.size(); ineuron++){
			layer[ineuron].CalculateInputError();
		}
	}

	// At this point, the input and output errors for all neurons in the ANN have been set.
	// We can now adjust the weights connecting the neurons using these values.

	// Loop over all connection layers adjusting the weights based on the errors
	for(unsigned int ilayer=0; ilayer<weights_layers.size(); ilayer++){
	
		// Loop overconnections in this connection layer and adjust the weights based
		// on the input to the connection (i.e. the output of the neuron in the preceeding layer)
		// and the error on the connection output (i.e. the input error of the neuron in the
		// following layer).
		nnann_weight_layer &weights = weights_layers[ilayer];
		for(unsigned int iweight=0; iweight<weights.size(); iweight++){
			nnann_weight &weight = weights[iweight];
			nnann_neuron *input = weight.GetInput();
			nnann_neuron *output = weight.GetOutput();

			// Standard backpropagation method for ANN adjusts the weight by an amount proportional
			// to the signal going into the connection and the error on the output end. That error
			// is, of course, cumulative containing contributions from all neurons in the previous
			// layer. However, it does make qualitative sense that the error on the weight would be
			// zero when the error on the thing it feeds is zero. It also makes sense that the weight
			// will have little influence if its input is zero. In that case, changing the weight won't
			// affect the input to the next layer so it's best to leave it alone. Making the error on the
			// weight propotional to these two is therefore not crazy. The learning rate "eta" limits how much
			// the weight changes for any given training sample requiring more samples/iterations (epochs)
			// to effect significant change. Thus, allowing the weights to represent the entire training
			// sample rather than just the last few patterns. (n.b. a similar mechanism exists for
			// Kalman filters.)
			double delta_weight = eta*input->GetOutput()*output->GetInputError();
//_DBG_<<"layer: "<<ilayer<<"  connection: "<<iweight<<"  weight:"<<weight.GetWeight()<<" delta_weight: "<<delta_weight<<" input's output:"<<input->GetOutput()<<" output's input(err):"<<output->GetInputError()<<" output's output(err):"<<output->GetOutputError()<<endl;
			weight.SetWeight(weight.GetWeight() + delta_weight);
		}
	}
	
	return sum_squared_error/(double)output_layer.size();
}

//---------------------------------
// FeedForward
//---------------------------------
void nnann::FeedForward(const std::vector<double> &inputs, std::vector<double> &outputs)
{
	// Dimension check
	if(layers.size()<1)throw nnann_exception("Empty network!");
	nnann_layer &input_layer = layers[0];

	if(inputs.size()!=input_layer.size())throw nnann_exception("Number of inputs don't match network!");

	// Clear the inputs for all neurons in the ANN
	for(unsigned int ilayer=0; ilayer<layers.size(); ilayer++)
		for(unsigned int ineuron=0; ineuron<layers[ilayer].size(); ineuron++)layers[ilayer][ineuron].ResetInput();

	// Set input layer to the input values
	for(unsigned int i=0; i<input_layer.size(); i++){
		// The first layer of neurons is just used to hold the inputs.
		// The activation function is not used for these
		input_layer[i].SetOutput(inputs[i]);
	}

	// Loop over all connection layers, propagating the signals forward
	for(unsigned int ilayer=0; ilayer<weights_layers.size(); ilayer++){
	
		// Loop over connections in this connection layer
		nnann_weight_layer &weights = weights_layers[ilayer];
		for(unsigned int iweight=0; iweight<weights.size(); iweight++){
			weights[iweight].ForwardPropagate();
		}
	
		// Loop over neurons in the layer we just propageted to and apply 
		// the activation function making their outputs valid
		vector<nnann_neuron> &layer = layers[ilayer+1];
		for(unsigned int ineuron=0; ineuron<layer.size(); ineuron++){
			layer[ineuron].CalculateOutput();
		}
	}
	
	// Copy outputs of network into provided container
	outputs.clear();
	nnann_layer &output_layer = layers[layers.size()-1];
	for(int ineuron=0; ineuron<output_layer.size(); ineuron++){
		outputs.push_back(output_layer[ineuron].GetOutput());
	}
}

//---------------------------------
// GetRMSError
//---------------------------------
double nnann::GetRMSError(std::vector<double> &target_vals)
{
	/// Get the RMS of the output layer errors based on the given target_vals.

	nnann_layer &output_layer = layers[layers.size()-1];
	double sum_diff2 = 0.0;
	for(unsigned int i=0; i<target_vals.size(); i++){
		double diff = target_vals[i] - output_layer[i].GetOutput();
		sum_diff2 += diff*diff;
	}

	return sqrt(sum_diff2/(double)target_vals.size());
}

//---------------------------------
// Initialize
//---------------------------------
void nnann::Initialize(void)
{
	/// This should be called if the network has not been initialized
	/// already by other means, prior to training. When called, all
	/// neurons for which an activation function is not set will
	/// be set to use a standard sigmoid activation function. In addition,
	/// all weights will be initialized to random values. 
	///
	/// This method is not called if the network has been loaded from
	/// another source (e.g. a file).

	// If activation function has not been set, then create a standard sigmoid function
	if(!activation_function)activation_function = new nnann_function_sigmoid();
	
	// Loop over all neruons, setting the activation function if not set
	for(unsigned int ilayer=1; ilayer<layers.size(); ilayer++){

		// Loop over neurons in the layer
		nnann_layer &layer = layers[ilayer];
		for(unsigned int ineuron=0; ineuron<layer.size(); ineuron++){
			nnann_neuron &neuron = layer[ineuron];
			
			// Set activation function if it's not already set
			if(!neuron.GetActivationFunction())neuron.SetActivationFunction(activation_function);			
		}
	}
	
	// Loop over all weights, setting them to random numbers
	for(unsigned int ilayer=0; ilayer<weights_layers.size(); ilayer++){
		nnann_weight_layer &weights = weights_layers[ilayer];
		for(unsigned int iweight=0; iweight<weights.size(); iweight++){
			double w = 2.0*(double)(random()&0x7FFFFFFF)/(double)(0x7FFFFFFF)-1.0;
			weights[iweight].SetWeight(w);
		}
	}

	initialized = true;
}

//---------------------------------
// LoadFromFile
//---------------------------------
void nnann::LoadFromFile(std::istream &ifs)
{
	// Delete any existing neurons from ann
	layers.clear();

	// If activation function has not been set, then create a standard sigmoid function
	if(!activation_function)activation_function = new nnann_function_sigmoid();

	// Read in layer info
	while(!ifs.eof() && !ifs.fail() && !ifs.bad()){
		string tok="";
		ifs>>tok;

		if(tok.length()>0 && tok[0]=='#'){
			// Comment line. Read in rest of line and continue
			char str[1024];
			ifs.getline(str,1024);
			continue;
		}else if(tok=="layer"){
			// Define a new layer (without weights)
			unsigned int ilayer, Nneurons;
			string function, colon;
			ifs >> ilayer >> colon >> Nneurons >> colon >> function;

			if(ilayer!=layers.size())throw nnann_exception("corrupt input file (layer number out of order)");
			
			// This will need to be changed to accomodate different activation functions
			nnann_function *my_function = NULL;
			if(function==activation_function->GetName())my_function=activation_function;
			
			// Create new layer with Nneurons and set the activation function for each
			AddLayer(Nneurons, my_function);

		}else if(tok=="Weights"){
			// Set the weights for an existing layer
			string forstr, layerstr;
			unsigned int ilayer;
			ifs >> forstr >> layerstr >> ilayer;
			
			if(ilayer<1 || ilayer>layers.size())throw nnann_exception("corrupt input file (layer number out of range)");

			nnann_weight_layer &weights = weights_layers[ilayer];
			for(unsigned int iweight = 0; iweight<weights.size(); iweight++){

				double w;
				ifs >> w;
				if(ifs.fail())throw nnann_exception("corrupt input file (unable to read weight!)");
				
				weights[iweight].SetWeight(w);
			}
			
			// If we just read in weights for the last layer, assume we are done. We should
			// stop reading from the input stream allowing any user extensions to be read
			// in. The file should be formatted such that the layer info is all read in first and
			// the weight info last.
			if(ilayer == (layers.size()-1))break;
		}
	}
	
	initialized = true;
}

//---------------------------------
// ReadFile
//---------------------------------
void nnann::ReadFile(std::string filename)
{
	/// Open a file with the specified file name and read the network in from it.
	/// Before closing the input file stream, call LoadExtensionsFromFile() so
	/// a sub-classes can read in additional, custom information at the end of the
	/// file.
	
	ifstream ifs(filename.c_str());
	if(!ifs.is_open()){
		throw nnann_exception(string("can't open \""+filename+"\"!"));
	}
	
	// Read core network from file
	LoadFromFile(ifs);

	// Call sub-class hook for reading in file.
	LoadExtensionsFromFile(ifs);
	
	ifs.close();
}

//---------------------------------
// RelinkLayers
//---------------------------------
void nnann::RelinkLayers(void)
{
	/// Create nnann_weight objects to link all of the neurons in the ANN.
	/// This will allocate the required nnann_weight objects to link all neurons
	/// between adjacent layers for the whole ANN. Input and output pointer lists 
	/// in the neurons will be updated to point to the new weight objects.
	///
	/// Note that this will preserve the actual weight values for nnann_weight objects
	/// that already exist (though the memory location of those objects may change). 
	/// New weights will be uninitialied however and so will need
	/// to either be set by hand or all the weights in the network set to random values
	/// by calling Initialize().

	// Less than 2 layers likely indicates an error
	if(layers.size()<2)throw nnann_exception("Less than 2 layers defined when nnann::Relink() called!");
	
	// Set the number of weights layers based on the content of layers
	weights_layers.resize(layers.size()-1);
	
	// Loop over pairs of layers starting with the first two
	for(unsigned int ilayer=0; ilayer<layers.size()-1; ilayer++){
		nnann_layer &prev = layers[ilayer];
		nnann_layer &next = layers[ilayer+1];
		nnann_weight_layer &weights = weights_layers[ilayer];
		
		// Resize the weights layer
		weights.resize(prev.size()*next.size() + next.size()); // extra next.size() is for connections to bias neuron
		
		// Clear the output weights list for all neurons in the previous layer
		for(unsigned int ineuron=0; ineuron<prev.size(); ineuron++) prev[ineuron].ClearOutputs();

		// Clear the input weights list for all neurons in the next layer
		for(unsigned int ineuron=0; ineuron<next.size(); ineuron++) next[ineuron].ClearInputs();

		// Loop over all pairs of neurons between the 2 layers and set
		// the pointers to the weight objects.
		for(unsigned int iprev=0; iprev<prev.size(); iprev++){
			nnann_neuron &prev_neuron = prev[iprev];
			for(unsigned int inext=0; inext<next.size(); inext++){
				nnann_neuron &next_neuron = next[inext];
				nnann_weight &weight = weights[iprev + inext*prev.size()];
				
				prev_neuron.AddOutput(&weight);
				next_neuron.AddInput(&weight);
				weight.SetInput(&prev_neuron);
				weight.SetOutput(&next_neuron);
			}
		}
		
		// Link the bias neuron to all "next" layer neurons
		for(unsigned int inext=0; inext<next.size(); inext++){
			nnann_neuron &next_neuron = next[inext];
			nnann_weight &weight = weights[prev.size()*next.size() + inext];
				
			bias.AddOutput(&weight);
			next_neuron.AddInput(&weight);
			weight.SetInput(&bias);
			weight.SetOutput(&next_neuron);
		}
	}
}

//---------------------------------
// SaveAs
//---------------------------------
void nnann::SaveAs(std::string filename)
{
	ofstream ofs(filename.c_str());
	if(!ofs.is_open()){
		throw nnann_exception(string("can't open \""+filename+"\"!"));
	}
	
	// Write header and core values to file
	WriteToFile(ofs);
	
	// Call sub-class hook for writing to file.
	WriteExtensionsToFile(ofs);
	
	ofs.close();
}

//---------------------------------
// TestOnFile
//---------------------------------
void nnann::TestOnFile(std::string filename)
{

}

//---------------------------------
// TrainFromFile
//---------------------------------
unsigned int nnann::TrainFromFile(std::string filename, double max_err, unsigned int max_epochs)
{
	// Make sure we're initialized
	if(!initialized)Initialize();
	
	// Make sure we have ANN of minimal topology
	if(layers.size()<2)throw nnann_exception("ANN has too few layers for training!");
	if(layers[0].size()<1)throw nnann_exception("ANN has too few inputs for training!");
	if(layers[layers.size()-1].size()<1)throw nnann_exception("ANN has too few outputs for training!");

	// Open input file
	ifstream ifs(filename.c_str());
	if(!ifs.is_open())throw nnann_exception(string("can't open \""+filename+"\"!"));
	
	// Read data from input file until we've read it all
	unsigned int Ninputs=0;
	unsigned int Noutputs=0;
	unsigned int epoch=0;
	while(!ifs.eof() && !ifs.fail() && !ifs.bad()){

		string tok="";
		ifs>>tok;

		if(tok.length()>0 && tok[0]=='#'){
			// Comment line. Read in rest of line and continue
			char str[1024];
			ifs.getline(str,1024);
			continue;
		}else if(tok=="inputs"){
			// Read in number of inputs
			ifs >> Ninputs;
		}else if(tok=="outputs"){
			// Read in number of outputs
			ifs >> Noutputs;
		}else if(tok=="data:"){
			// Check that Ninputs and Noutputs are are correct
			if(Ninputs!=layers[0].size())throw nnann_exception(string("Number of inputs in training file doesn't match ANN!"));
			if(Noutputs!=layers[layers.size()-1].size())throw nnann_exception(string("Number of outputs in training file doesn't match ANN!"));
			
			// Remember stream position pointer so we can loop
			streampos pos_start = ifs.tellg();
			
			cout<<"Training on \""<<filename<<"\" with "<<Ninputs<<" inputs and "<<Noutputs<<" ouputs ..."<<endl;
			
			// Loop over entire training set repeatedly until we either reach max_epochs
			// or achieve the desired accuracy.
			for(epoch=0; epoch<max_epochs||max_epochs==0; epoch++){
				ifs.clear();
				ifs.seekg(pos_start);
				double err = TrainOnePass(ifs);
				if(err<max_err)break;
				if((epoch%1000) == 0)cout<<"err="<<err<<" after "<<epoch<<" epochs"<<endl;
			}
		}
	}
	
	// Close input file
	ifs.close();
	
	return epoch;
}

//---------------------------------
// TrainOnePass
//---------------------------------
double nnann::TrainOnePass(ifstream &ifs)
{
	unsigned int Ninputs = layers[0].size();
	unsigned int Noutputs = layers[layers.size()-1].size();

	// Loop over training samples
	double sum_diff2 = 0.0;
	double Ndiff2 = 0.0;
	while(true){
		vector<double> inputs(Ninputs);
		vector<double> target_vals(Noutputs);
		for(unsigned int i=0; i<Ninputs; i++)ifs >> inputs[i];
		for(unsigned int i=0; i<Noutputs; i++)ifs >> target_vals[i];
		if(ifs.fail() || ifs.eof())break;
		
		// Run the network using the input values
		vector<double> outputs;
		FeedForward(inputs, outputs);
		
		// Backpropagate based on target values, updating weights.
		double diff2 = BackPropagate(target_vals);
		sum_diff2 += diff2;
		Ndiff2+=1.0;

#if 0
		cout<<"Training inputs:";
		for(unsigned int i=0; i<Ninputs; i++)cout<<" "<<inputs[i];
		cout<<endl;
		cout<<"Training outputs:";
		for(unsigned int i=0; i<Noutputs; i++)cout<<" "<<target_vals[i];
		cout<<endl;
		cout<<"     ANN outputs:";
		for(unsigned int i=0; i<Noutputs; i++)cout<<" "<<outputs[i];
		cout<<endl;
		cout<<"RMS: "<<sqrt(diff2)<<endl;
		cout<<endl;
#endif
	}
	
	if(Ndiff2==0.0)throw nnann_exception("Empty training set!");
	
	double err = sqrt(sum_diff2/Ndiff2);

#if 0
	cout<<"============ Total error: "<<err<<" =================="<<endl;
#endif

	return err;
}

//---------------------------------
// WriteToFile
//---------------------------------
void nnann::WriteToFile(std::ostream &ofs)
{
	// Write header info
	ofs<<"#"<<endl;
	ofs<<"#"<<endl;
	ofs<<"#"<<endl;
	ofs<<endl;
	
	// Write number of neurons per layer and activation function used
	// by neurons in the layer.
	ofs<<"# layer : Num. neurons : activation function"<<endl;
	for(unsigned int ilayer=0; ilayer<layers.size(); ilayer++){
		nnann_function *activation_function = layers[ilayer][0].GetActivationFunction();
		string function_name = activation_function ? activation_function->GetName():"Unknown";
		ofs<<"layer "<<ilayer<<" : "<<layers[ilayer].size()<<" : "<<function_name<<endl;
	}
	ofs<<endl;
	ofs<<endl;
	
	// Write weights for each layer
	for(unsigned int ilayer=0; ilayer<weights_layers.size(); ilayer++){
		ofs<<"Weights for layer "<<ilayer<<endl;
		
		nnann_weight_layer &weights = weights_layers[ilayer];
		unsigned int Nneurons = layers[ilayer].size();
		for(unsigned int iweight=0; iweight<weights.size(); iweight++){
			
			ofs << weights[iweight].GetWeight() << " ";
			if((iweight%Nneurons) == (Nneurons-1))ofs<<endl;

		}
		ofs<<endl;
	}	
}