/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2002 Markus Turtinen and Topi Mäenpää
 * All rights reserved.
 *
 * This program is free software. You can redistribute and/or modify
 * it under the terms of the free software licence found in the
 * accompanying file "COPYING". The licence terms must always be
 * redistributed with this source file. The above copyright notice
 * must be reproduced in all modified and unmodified copies of this
 * source file.
 *
 * $Revision: 1.6 $
 *********************************************************************/


#ifndef _SOM_H
#define _SOM_H

#include "../VectorQuantizer.h"
#include "../Random.h"
#include <ListUtils.h>
#include <io/IO.h>
#include <fstream>

namespace prapi { namespace neuro {

	/**
	 * SOM topologies.
	 *
	 * <ul>
	 * <li>SOM_HEXAGONAL - six neighbors have a distance of one</li>
	 * <li>SOM_SQUARE - four neighbors have a distance of one</li>
	 * </ul>
	 **/
	enum SOMTopology { SOM_HEXAGONAL, SOM_SQUARE };

	/**
	 * SOM learning rate functions.
	 *
	 * <ul>
	 * <li>SOM_LINEAR_ALPHA - learning rate decreases linearly</li>
	 * <li>SOM_INVERSE_T_ALPHA - learning rate is inversely proportional to training interation index</li>
	 * </ul>
	 **/
	enum SOMLearningRate { SOM_LINEAR_ALPHA, SOM_INVERSE_T_ALPHA };

	/**
	 * A base class for different types of neighborhoods. When updating
	 * nodes in the SOM, the amount of vector movement is determined by
	 * the weight returned by a SOMNeighborhood instance and the current
	 * learning rate (multiplied).
	 **/
	class SOMNeighborhood : virtual public util::Object
	{
	public:
		/**
		 * Get the weight for a node update.
		 *
		 * @param radius the current effective update radius
		 * @param distance the distance from the winning node to the one
		 * to be updated.
		 * @param the weight for a vector update, [0,1]. 0 means no change.
		 **/
		virtual double getWeight(double radius, double distance) const = 0;

		class Bubble;
		class Gaussian;
	};

	/**
	 * "Bubble" neighborhood. Each node within the current radius is
	 * updated with a weight of one. Others are not updated.
	 **/
	class SOMNeighborhood::Bubble : public SOMNeighborhood
	{
	public:
		double getWeight(double radius, double distance) const { return distance <= radius ? 1 : 0; }
	};

	/**
	 * An implementation of the self-organizing map (Kohonen map). 
	 * Typically, SOM is trained in two phases. First, a "rough"
	 * training is made with a relatively large learning rate (~0.05), a
	 * relatively large radius (e.g. 10, depending on the size of the
	 * map), and a relatively low number of iterations (depending on the
	 * number of training samples). Then, fine-tuning is performed with
	 * a smaller learning rate (~0.02), a smaller radius (e.g. 3), and a
	 * larger number of iterations.
	 **/
	template <class I=std::string, class C=int> class SOM : public VQClassifier<I,C>
	{
	public:
		/**
		 * Create a new SOM with the given size.
		 **/
		SOM(int xSize=1, int ySize=1) :
			VQClassifier<I,C>(new EuclideanDistance<double>),
			_iSizeX(xSize), _iSizeY(ySize),
			_lLearningLength(1000), _dRadius(10), _dLearningRate(0.05),
			_topology(SOM_HEXAGONAL),	_learningRateFunction(SOM_LINEAR_ALPHA),
			_pNeighborhood(new SOMNeighborhood::Bubble) {}


		/**
		 * Initialize the code book. The number of code vectors must be
		 * equal to the size of the SOM, i.e. width*height. Typically, the
		 * initial code book consists of samples selected randomly from
		 * training data.
		 **/
		void initialize(const util::List<Sample<double,I,C> >& initialCodebook)
			throw (InvalidArgumentException&);

		/**
		 * Train the SOM with the given training data.
		 **/
		void train(const util::List<Sample<double,I,C> >& trainingSamples)
			throw (ClassificationException&);
		
		/**
		 * Give class labels to SOM nodes. A simple voting procedure is
		 * used here: each training sample is assigned to the closest
		 * node, and each node gets a class label of the most common class
		 * in it.
		 **/
		void label(const util::List<Sample<double,I,C> >& trainingSamples)
			throw (ClassificationException&);

		/**
		 * Set the learning radius. The radius affects the number of code
		 * vectors that are adapted in the neighborhood of the winning
		 * node.
		 **/
		void setRadius(double radius) { _dRadius = radius; }
		/**
		 * Get the learning radius.
		 **/
		double getRadius() const { return _dRadius; }
		/**
		 * Set the topology of the network.
		 **/
		void setTopology(SOMTopology topology) { _topology = topology; }
		/**
		 * Get the topology of the network.
		 **/
		SOMTopology getTopology() const { return _topology; }
		/**
		 * Set the number of iterations in training.
		 **/
		void setLearningLength(long length) { _lLearningLength = length; }
		/**
		 * Set learning rate.
		 **/
		void setLearningRate(double rate) { _dLearningRate = rate; }
		/**
		 * Set the type of learning rate change in training.
		 **/
		void setLearningRateFunction(SOMLearningRate func) { _learningRateFunction = func; }
		/**
		 * Get the number of iterations in training.
		 **/
		long getLearningLength() const { return _lLearningLength; }
		/**
		 * Get learning rate.
		 **/
		double getLearningRate() const { return _dLearningRate; }
		/**
		 * Get the type of learning rate change in training.
		 **/
		double getLearningRateFunction() const { return _learningRateFunction; }
		/**
		 * Get the width of the SOM.
		 **/
		int getWidth() const { return _iSizeX; }
		/**
		 * Get the height of the SOM.
		 **/
		int getHeight() const { return _iSizeY; }

		virtual ~SOM() {}

		/**
		 * Write a SOM_PAK formatted SOM description to a file. SOM_PAK is
		 * a well-known Matlab implementation of the SOM.
		 *
		 * @param fileName the name of the file to write into
		 * @param neighborhoodName the type of the neighborhood function
		 * represented textually. Use only if you know what you are doing.
		 * @param classNames a list of textual names for the class
		 * indices. If this list is empty, class indices will be used as
		 * class names.
		 * @exception IOException if the file cannot be written to
		 **/
		void writeSOM_PAKFile(std::string fileName,
													std::string neighborhoodName = "bubble",
													const util::List<std::string>& classNames = util::List<std::string>(0))
			throw (util::io::IOException&);

	private:
		int _iSizeX, _iSizeY;
		long _lLearningLength;
		double _dRadius, _dLearningRate;
		SOMTopology _topology;
		SOMLearningRate _learningRateFunction;
		util::SmartPtr<SOMNeighborhood> _pNeighborhood;
		
		double linearAlpha(long iter)
		{
			return (_dLearningRate * (double) (_lLearningLength - iter) / (double) _lLearningLength);
		}
		
		double inverseTAlpha(long iter)
		{
			double c = (double)_lLearningLength / 100.0;
			return (_dLearningRate * c / (c + iter));
		}
		
		double distSquare(int bx, int by, int tx, int ty);
		double distHexa(int bx, int by, int tx, int ty);

		void adaptNodes(int hitX, int hitY, const Sample<double,I,C>& sample,
										double radius, double alpha);

	protected:
		/**
		 * Adapt a code vector towards a sample. The default
		 * implementation calculates a weighed average of the sample and
		 * code vectors.
		 *
		 * @param code the code vector to be adapted
		 * @param sample the training sample according to which the code
		 * vector is to be adapted
		 * @param alpha adaptation factor, obtained by multiplying the
		 * current learning rate by the value of the neighborhood function
		 * at the position of the code vector.
		 **/
		virtual void adaptVector(Sample<double,I,C>& code, const Sample<double,I,C>& sample, double alpha);
	};

	
	template <class I, class C> void SOM<I,C>::train(const util::List<Sample<double,I,C> >& trainingSamples)
		throw (ClassificationException&)
	{
		//Rearrange samples randomly
		int len = trainingSamples.getLength();
		List<int> randomIndices(len);
		for (int i=0; i<len; i++)
			randomIndices += i;
		Random::init();
		Random::shuffle(randomIndices);
		
		//Training loop
    for(long int i=0; i<_lLearningLength; i++)
			{
				//radius decreases linearly from initial value to 1
        double nextRad = (double)(1.0 + (_dRadius-1.0) * (double)(_lLearningLength-i)/(double)_lLearningLength);
				double nextAlpha = 0;

        //teaching rate decreases linearly from initial value to 0
        if(_learningRateFunction == SOM_LINEAR_ALPHA)
					nextAlpha = linearAlpha(i);
        else 
					nextAlpha = inverseTAlpha(i);
				
        int hitnum = getBinIndex(trainingSamples[randomIndices[i%len]]);
        int hX = hitnum % _iSizeX;
        int hY = hitnum / _iSizeX;

				adaptNodes(hX, hY, trainingSamples[randomIndices[i%len]], nextRad, nextAlpha);
			}
	}

	template <class I, class C> void SOM<I,C>::adaptNodes(int hitX, int hitY, const Sample<double,I,C>& sample,
																												double radius, double alpha)
	{
    for(int index=0; index<_lstCodeBook.getLength(); index++)
			{
				//current x,y
        int tX = index % _iSizeX;
        int tY = index / _iSizeX;

				double weight;
        if(_topology = SOM_HEXAGONAL)
					weight = _pNeighborhood->getWeight(radius, distHexa(hitX, hitY, tX, tY));
				else
					weight = _pNeighborhood->getWeight(radius, distSquare(hitX, hitY, tX, tY));

				if (weight != 0)
					adaptVector(_lstCodeBook[index], sample, alpha*weight);
			}
	}

	//distance between two nodes in hexagonal topology
	template <class I, class C> double SOM<I,C>::distHexa(int bx, int by, int tx, int ty)
	{
		double ret, diff;
    diff = bx - tx;
    if (((by - ty) % 2) != 0)
			{
				if ((by % 2) == 0)
					diff -= 0.5;
				else
					diff += 0.5;
			}
    ret = diff * diff;
    diff = by - ty;
    ret += 0.75 * diff * diff;
    return(sqrt(ret));
	}

	//distance between two nodes in rectangular topology
	template <class I, class C> double SOM<I,C>::distSquare(int bx, int by, int tx, int ty)
	{
		double ret, diff;
    diff = bx - tx;
    ret = diff * diff;
    diff = by - ty;
    ret += diff * diff;
    return(sqrt(ret));
	}

	template <class I, class C> void SOM<I,C>::adaptVector(Sample<double,I,C>& code, const Sample<double,I,C>& sample, double alpha)
	{
		double tmp = 1.0-alpha;
		for (int i=code.featureVector().getLength(); i--;)
			code.featureVector()[i] = alpha * sample.featureVector()[i] + tmp * code.featureVector()[i];
	}

	template <class I, class C> void SOM<I,C>::initialize(const util::List<Sample<double,I,C> >& codeBook)
		throw (InvalidArgumentException&)
	{
		if (codeBook.getLength() != _iSizeX * _iSizeY)
			throw InvalidArgumentException("SOM::initialize(List<Sample>&): Must have " +
																		 String::toString(_iSizeX*_iSizeY) + " code vectors.");
		_lstCodeBook = codeBook;
	}

	template <class I, class C> void SOM<I,C>::label(const util::List<Sample<double,I,C> >& trainingSamples)
		throw (ClassificationException&)
	{
		util::List<util::List<int> > hits(_lstCodeBook.getLength());
		hits.setLength(_lstCodeBook.getLength());
		int classCount = 0;

		//Go through training samples and place class index votes on
		//winning nodes
		for(int i=0; i<trainingSamples.getLength(); i++)
			{
				int bestNode = getBinIndex(trainingSamples[i]);
				int classIndex = (int)trainingSamples[i].getTrueClass();
				if (hits[bestNode].getLength() <= classIndex)
					hits[bestNode].setLength(classIndex+1, 0);
				hits[bestNode][classIndex]++;
				if (classIndex > classCount)
					classCount = classIndex;
			}
		classCount++;

		//Find winner classes
		for (int i=_lstCodeBook.getLength(); i--;)
			{
				int trueClass = -1;
				if (hits[i].getLength())
					trueClass = ListUtils::maxIndex(hits[i]);
				_lstCodeBook[i].trueClass() = trueClass;
			}

		//Eliminate empty nodes
		kNNClassifier<double,I,C> knn(const_cast<util::List<Sample<double,I,C> >&>(trainingSamples), *_proximityMeasure, classCount, 3);
		for (int i=_lstCodeBook.getLength(); i--;)
			if ((int)_lstCodeBook[i].trueClass() == -1)
				_lstCodeBook[i].trueClass() = knn.getClassification(_lstCodeBook[i]);
	}

	template <class I, class C>	void SOM<I,C>::writeSOM_PAKFile(std::string fileName,
																															std::string neighborhoodName,
																															const util::List<std::string>& classNames)
		throw (util::io::IOException&)
	{
		std::ofstream out(fileName.c_str());
		if (!out)
			throw util::io::IOException("SOM::writeSOM_PAKFile(string, string): cannot open " + fileName + ".");
		out << _lstCodeBook[0].featureVector().getLength() << " " << (_topology == SOM_HEXAGONAL ? "hexa " : "rect ")
				<< _iSizeX << " " << _iSizeY << " " << neighborhoodName << endl;
		for (int i=0; i<_lstCodeBook.getLength(); i++)
			{
				for (int j=0; j<_lstCodeBook[i].featureVector().getLength(); j++)
					{
						if (j) out << " ";
						out << _lstCodeBook[i].featureVector()[j];
					}

				if (_lstCodeBook[i].trueClass() >= 0)
					{
						out << " ";
						if (classNames.getLength())
							out << classNames[_lstCodeBook[i].trueClass()] << endl;
						else
							out << _lstCodeBook[i].trueClass() << endl;
					}
			}
		out.close();
	}
}}

#endif