/*********************************************************************
 * This file is part of the cpplibs suite.
 *
 * Copyright (C) 2001 Topi Mäenpää and Jaakko Viertola
 * 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.24 $
 *********************************************************************/

#ifndef _MATHCLASS_H
#define _MATHCLASS_H

#include <math.h>
#include "Matrix.h"
#include "Heap.h"
#include "Exception.h"

namespace util
{
	/**
	 * MathException is thrown in many cases where invalid mathematic
	 * operations would be performed.
	 **/
	class MathException : public Exception
	{
	public:
		MathException(std::string message) : Exception(message) {}
	};
	
	/**
	 * A template class that holds some useful functions for calculating
	 * mathematic quantities. Note that these operations are supposed to
	 * work with elementary types only.
	 **/
	class Math
	{
	public:
		/**
		 * Possible types of variance. When calculating variance in
		 * VAR_BIASED mode, the result is
		 * 1/N*sum<sub>i=1..N</sub>(x<sub>i</sub>-u<sub>x</sub>)<sup>2</sup>.
		 * In VAR_UNBIASED mode the divisor 1/N is changed to 1/(N-1) to
		 * obtain an unbiased estimate of the sample variance. Stddev uses
		 * VAR_UNBIASED by default.
		 **/
		enum VarianceType { VAR_BIASED=0, VAR_UNBIASED=1 };
		/**
		 * Calculate the variance of all values in a list. The result
		 * value is the best unbiased estimate of var(x), i.e. 1/(n-1) *
		 * sum((x-m<sub>x</sub>)^2).
		 **/
		template <class T> static double variance(const List<T>& lst, VarianceType type = VAR_UNBIASED) throw (MathException&);
		/**
		 * Calculate the standard deviation of all values in a list.
		 **/
		template <class T> static double stdev(const List<T>& lst) throw (MathException&);
		/**
		 * Calculate the mean of all values in a list.
		 **/
		template <class T> static double mean(const List<T>& lst) throw (MathException&);
		/**
		 * Find the median in a list of values.
		 **/
		template <class T> static T median(const List<T>& lst) throw (MathException&);
		/**
		 * Calculate the mean and variance for the values in a list.
		 **/
		template <class T> static void meanAndVariance(const List<T>& mat,
																									 double& mean,
																									 double& variance,
																									 VarianceType type = VAR_UNBIASED) throw (MathException&);
		/**
		 * Finds the min and max values from a matrix.
		 **/ 
		template <class T> static void minAndMax(const Matrix<T>& lst, T& min, T& max) throw (MathException&);

		/**
		 * Calculate the factorial of x (i.e. x!). For large values of x
		 * the returned value may not be quite accurate.
		 **/
		static double factorial(unsigned int x);
		/**
		 * Calculate the sum of all values in a list.
		 **/
		template <class T> static T sum(const List<T>& lst) throw (MathException&);
		/**
		 * Finds the maximum value of the list.
		 **/
		template <class T> static T max(const List<T>& lst) throw (MathException&);
		/**
		 * Finds the minimum value of the list.
		 **/
		template <class T> static T min(const List<T>& lst) throw (MathException&);
		/**
		 * Finds the min and max values from a list.
		 **/ 
		template <class T> static void minAndMax(const List<T>& lst, T& min, T& max) throw (MathException&);
		/**
		 * Find the number of given parameters from list.
		 **/
		template <class T> static int numberOf(const List<T>& lst,const T parameter) throw (MathException&);
		/**
		 * Calculate the entropy of the elements in a list.
		 **/
		template <class T> static double entropy(const List<T>& lst) throw (MathException&);
		
		/**
		 * Calculate the variance of all values in a matrix.
		 **/
		template <class T> static double variance(const Matrix<T>& mat, VarianceType type = VAR_UNBIASED) throw (MathException&);
		/**
		 * Calculate a covariance matrix. In the input matrix, each row is
		 * an observation and each column a variable. The diagonal of the
		 * covariance matrix contains variances for each variable. The
		 * resulting matrix will be of size NxN where N is the number of
		 * variables (dimensions) in each observation vector.
		 **/
		template <class T> static Matrix<double> covariance(const Matrix<T>& mat) throw (MathException&,MatrixException&);
		/**
		 * Calculate the standard deviation of all values in a matrix.
		 **/
		template <class T> static double stdev(const Matrix<T>& mat) throw (MathException&);
		/**
		 * Calculate the mean of all values in a matrix.
		 **/
		template <class T> static double mean(const Matrix<T>& mat) throw (MathException&);
		/**
		 * Find the median in a matrix.
		 **/
		template <class T> static T median(const Matrix<T>& mat) throw (MathException&);
		/**
		 * Calculate mean and variance for all entries in a matrix.
		 **/
		template <class T> static void meanAndVariance(const Matrix<T>& mat,
																									 double& mean,
																									 double& variance,
																									 VarianceType type = VAR_UNBIASED) throw (MathException&);
		/**
		 * Calculate the sum of all values in a matrix.
		 **/
		template <class T> static T sum(const Matrix<T>& mat) throw (MathException&) { return mat.sum(); }
		/**
		 * Find the maximum value in a matrix.
		 **/
		template <class T> static T max(const Matrix<T>& mat) throw (MathException&) { return mat.max(); }
		/**
		 * Find the minimum value in a matrix.
		 **/
		template <class T> static T min(const Matrix<T>& mat) throw (MathException&) { return mat.min(); }
		/**
		 * Find the number of given parameters from the matrix.
		 **/
		template <class T> static int numberOf(const Matrix<T>& mat,const T& parameter) throw (MathException&);
		
		/**
		 * Calculate the entropy of the elements in a matrix.
		 **/
		template <class T> static double entropy(const Matrix<T>& mat) throw (MathException&);
		/**
		 * Calculate the combination given in whichCombinations
		 * from classes given in numberOfClasses (whis is the MAX index
		 * for class. fl. if you have classes 0,1,2 give number 2 and
		 * function calculates combinations).
		 * @param numberOfClasses the max index of classes
		 * @param whichCombinations the combination wanted
		 **/
		static List<List<int> > findCombinations(int numberOfClasses,int whichCombinations)
			throw (MathException&);
		/**
		 * Calculate all combinations of given number of classes.
		 * @param numberOfClasses the max index of classes
		 **/
		static List<List<int> > findAllCombinations(int numberOfClasses)
			throw (MathException&);
	};

	template <class T> double Math::entropy(const List<T>& lst) throw (MathException&)
	{
		double s = sum(lst);
		double entropy = 0;
		for (int i=lst.getLength();i--;)
			{
				double tmp = double(lst[i])/s;
				entropy -= tmp*log(tmp);
			}
		return entropy;
	}
	
	template <class T> double Math::variance(const List<T>& lst, VarianceType type) throw (MathException&)
	{
		const T* data = lst.getData();
		int len = lst.getLength();
		double sum = 0, sSum = 0;
		for (int i=0;i<len;i++,data++)
			{
				sum += double(*data);
				sSum += double(*data)*double(*data);
			}
		return (sSum - sum*sum/len) / (len-(int)type);
	}

	template <class T> void Math::meanAndVariance(const List<T>& lst,
																								double& mean,
																								double& variance,
																								VarianceType type) throw (MathException&)
	{
		const T* data = lst.getData();
		int len = lst.getLength();
		double sum = 0, sSum = 0;
		for (int i=0;i<len;i++,data++)
			{
				sum += double(*data);
				sSum += double(*data)*double(*data);
			}
		mean = sum/len;
		variance = (sSum - sum*mean) / (len-(int)type);
	}

	template <class T> double Math::stdev(const List<T>& lst) throw (MathException&)
	{
		return sqrt(variance(lst));
	}
	
	template <class T> double Math::mean(const List<T>& lst) throw (MathException&)
	{
		int len = lst.getLength();
		if (len > 0)
			{
				double sum = 0;
				for (int i=0;i<len;i++)
					sum += (double)lst[i];
				return sum/len;
			}
		else
			throw MathException("Math::mean(const List<T>&): Cannot calculate mean for an empty list.");
	}

	template <class T> T Math::median(const List<T>& lst) throw (MathException&)
	{
		int len = lst.getLength();
		if (len > 0)
			{
				Heap<T> heap(len);
				heap.addElements(lst);
				for (int i=0;i<len>>1;i++)
					heap.removeElementAt(0);
				return heap[0];
			}
		else
			throw MathException("Math::median(const List<T>&): Cannot find median in an empty list.");
	}

	template <class T> T Math::sum(const List<T>& lst) throw (MathException&)
	{
		int len = lst.getLength();
		if (len > 0)
			{
				T sum = lst[0];
				for(int i=1;i<len;i++)sum += lst[i];
				return sum;
			}
		else
			throw MathException("Math::sum(const List<T>&): Cannot calculate sum for an empty list.");
	}

	template <class T> T Math::max(const List<T>& lst) throw (MathException&)
	{
		int len = lst.getLength();
		if (len > 0)
			{
				T max = lst[0];
				for(int i=1;i<len;i++)
					if(max < lst[i])max = lst[i];
				return max;
			}
		else
			throw MathException("Math::max(const List<T>&): Cannot calculate max value for an empty list.");
	}
		
	template <class T> T Math::min(const List<T>& lst) throw (MathException&)
	{
		int len = lst.getLength();
		if (len > 0)
			{
				T min = lst[0];
				for(int i=1;i<len;i++)
					if(lst[i] < min)min = lst[i];
				return min;
			}
		else
			throw MathException("Math::min(const List<T>&): Cannot calculate minium value for an empty list.");
	}
	
	template <class T> void Math::minAndMax(const List<T>& lst, T& min, T& max) throw (MathException&)
	{
		int len = lst.getLength();
		if (len > 0)
			{
				min = lst[0];
				max = lst[0];
				for(int i=1;i<len;i++)
					{
						if(lst[i] < min)min = lst[i];
						else if(max < lst[i])max = lst[i];
					}
			}
		else
			throw MathException("Math::minAndMax(const List<T>&): Cannot calculate minium and maximum value for an empty list.");		
	}

	template <class T> void Math::minAndMax(const Matrix<T>& mat, T& min, T& max) throw (MathException&)
	{
		const T* data = mat.getData();
		int size = mat.getRows() * mat.getColumns();
		if (size > 0)
			{
				min = *data;
				max = *data++;
				for (int i=size-1; i--; data++)
					{
						if(*data < min) min = *data;
						else if(*data > max) max = *data;
					}
			}
		else
			throw MathException("Math::minAndMax(const Matrix<T>&): Cannot calculate minium and maximum value for an empty matrix.");		
	}

	template <class T> int Math::numberOf(const List<T>& lst,const T parameter) throw (MathException&)
	{
		int len = lst.getLength();
		int result =0;
		if (len > 0)
			{
				for(int i=0;i<len;i++)
					if(lst[i] ==  parameter)result++;
				return result;
			}
		else
			throw MathException("Math::numberOf(const List<T>&): Cannot calculate numberOf value for an empty list.");	
	}

	template <class T> void Math::meanAndVariance(const Matrix<T>& mat,
																								double& mean,
																								double& variance,
																								VarianceType type) throw (MathException&)
	{
		const T* data = mat.getData();
		int size = mat.getRows() * mat.getColumns();
		double sum = 0, sSum = 0;
		for (int i=0;i<size;i++,data++)
			{
				sum += double(*data);
				sSum += double(*data)*double(*data);
			}
		mean = sum/size;
		variance = (sSum - sum*mean) / (size-(int)type);
	}
	
	template <class T> double Math::entropy(const Matrix<T>& mat) throw (MathException&)
	{
		const T* data = mat.getData();
		int size = mat.getRows() * mat.getColumns();
		double s = mat.sum();
		double entropy = 0;
		for (int i=size;i--;data++)
			{
				double tmp = double(*data)/s;
				entropy -= tmp*log(tmp);
			}
		return entropy;
	}

	template <class T> double Math::variance(const Matrix<T>& mat, VarianceType type) throw (MathException&)
	{
		const T* data = mat.getData();
		int size = mat.getRows() * mat.getColumns();
		double sum = 0, sSum = 0;
		for (int i=0;i<size;i++,data++)
			{
				sum += double(*data);
				sSum += double(*data)*double(*data);
			}
		return (sSum - sum*sum/size) / (size-(int)type);
	}

	template <class T> Matrix<double> Math::covariance(const Matrix<T>& mat) throw (MathException&,MatrixException&)
	{
		int dims = mat.getColumns();
		int samples = mat.getRows();
		List<double> means(dims);
		means.setLength(dims,0);
		for (int i=0;i<dims;i++)
			{
				for (int j=0;j<samples;j++)
					means[i] += mat(j,i);
				means[i] /= samples;
			}

		const T* matData = mat.getData();
		Matrix<double> result(dims,dims);
		Matrix<double> diff(dims,1);
		for (int i=0;i<samples;i++)
			{
				double* diffData = diff.getData();
				for (int j=0;j<dims;j++,diffData++,matData++)
					*diffData = *matData - means[j];
				result += diff * diff.getTranspose();
			}
		result /= (samples-1);
		return result;
	}
	
	template <class T> double Math::stdev(const Matrix<T>& mat) throw (MathException&)
	{
		return sqrt(variance(mat));
	}
	
	template <class T> double Math::mean(const Matrix<T>& mat) throw (MathException&)
	{
		int size = mat.getRows() * mat.getColumns();
		if (size > 0)
			return (double)mat.sum()/(double)size;
		else
			throw MathException("Math::mean(const Matrix<T>&): Cannot calculate mean for an empty matrix.");
	}
	
	template <class T> T Math::median(const Matrix<T>& mat) throw (MathException&)
	{
		int size = mat.getRows() * mat.getColumns();
		if (size > 0)
			{
				Heap<T> heap(size);
				const T* ptr = mat.getData();
				for (int i=size;i--;ptr++)
					heap += *ptr;
				
				for (int i=0;i<size>>1;i++)
					heap.removeElementAt(0);
				return heap[0];
			}
		else
			throw MathException("Math::median(const Matrix<T>&): Cannot find median in an empty matrix.");
	}

	template <class T> int Math::numberOf(const Matrix<T>& mat,const T& parameter) throw (MathException&)
	{
		int len = mat.getRows()*mat.getColumns();
		const T* data = mat.getData();
		int result =0;
		if (len > 0)
			{
				for(int i=0;i<len;i++,data++)
					if(*data ==  parameter)result++;
				return result;
			}
		else
			throw MathException("Math::numberOf(const Matrix<T>&): Cannot calculate numberOf value for an empty Matrix.");	
	}
}

#endif