/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2001-2002 Topi Mäenpää
 * Copyright (C) 2001 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.25 $
 *********************************************************************/

#ifndef _PROXIMITYMEASURE_H
#define _PROXIMITYMEASURE_H

#include <math.h>
#include <values.h>
#include <List.h>
#include <Matrix.h>
#include <Math.h>
#include <string>
#include "Sample.h"
#include "Cluster.h"
#include "ClassificationException.h"

namespace prapi
{
	/**
	 * ProximityException is used if proximity calculation cannot be
	 * performed for some reason.
	 **/
	class ProximityException : public ClassificationException
	{
	public:
		ProximityException(std::string msg) : ClassificationException(msg) {}
	};

	/**
	 * ProximityMeasure is a general representation of a proximity
	 * measure between samples, sample sets (clusters) or between a
	 * sample and a cluster. It is used by classifiers and the
	 * clustering algorighms provided by the clustering library. In
	 * these contextes, "proximity" is always taken to mean
	 * "dissimilarity". That is, minimizing the proximity measure means
	 * finding closest sample or cluster. If you are using similarity
	 * measure, try to find a reasonable inversion for it.
	 **/
	template <class T> class ProximityMeasure : virtual public util::Object
	{
	public:
		/**
		 * Create a new ProximityMeasure that either uses or does not use
		 * cluster representatives and either is or is not symmetric.
		 **/
		ProximityMeasure(bool representatives=true, bool symmetric=true) :
			_bUseRepresentatives(representatives), _bSymmetric(symmetric) {}
		virtual ~ProximityMeasure() {}

		/**
		 * Proximity between two lists.
		 *
		 * @param sample the list to be classified
		 * @param model the model
		 * @param stopAfter A value that can be used in optimizing the
		 * performance of a proximity measure. Classifiers may set this
		 * value to inform the proximity measure to stop if it notices the
		 * proximity will be over a certain value. For example, the kNN
		 * classifier sets this value to the <i>k</i>th smallest proximity
		 * so far. If a sample is not closer than that, then it makes no
		 * sense to inspect it any more.
		 **/
		virtual double getProximity(const util::List<T>& sample, const util::List<T>& model,
																double stopAfter = MAXDOUBLE) const
			throw (ProximityException&) = 0;

		/**
		 * Proximity between two samples.
		 * @param sample the sample to be classified
		 * @param model the model
		 * @param stopAfter maximum allowable proximity
		 * @see #getProximity(const util::List&, const util::List&, double)
		 **/
		template <class I,class C> double getProximity(const Sample<T,I,C>& sample,
																									 const Sample<T,I,C>& model,
																									 double stopAfter = MAXDOUBLE) const
			throw (ProximityException&);
		/**
		 * Proximity between sample and a set. If cluster representatives
		 * are set off, this implementation returns the minimum proximity
		 * between the given sample and any sample in the given cluster by
		 * using the getProximity method. Otherwise the proximity between
		 * the sample and the cluster's representative is returned.
		 *
		 * If representatives are not used, the complexity of this
		 * procedure is O(N), where N is the number of samples in a
		 * cluster.
		 *
		 * @param sample the sample to be inspected
		 * @param sampleSet a set of samples that form a cluster
		 * @param stopAfter maximum allowable proximity
		 * @see #getProximity(const util::List&, const util::List&, double)
		 **/
 		template <class I,class C> double getProximity(const Sample<T,I,C>& sample,
																									 const Cluster<T,I,C>& sampleSet,
																									 double stopAfter = MAXDOUBLE) const
			throw (ProximityException&);
		/**
		 * Proximity between two sets. If cluster representatives are not
		 * used, this implementation returns the minimum proximity between
		 * any two samples in the given sets by using the getProximity
		 * method. Otherwise the proximity between cluster representatives
		 * is returned.
		 *
		 * If representatives are not used, the complexity of this
		 * procedure is O(NM), where N and M are the sample counts in each
		 * cluster
		 *
		 * @param sampleSet1 cluster one
		 * @param sampleSet1 cluster two
		 * @param stopAfter maximum allowable proximity
		 * @see #getProximity(const util::List&, const util::List&, double)
		 **/
	 	template <class I,class C> double getProximity(const Cluster<T,I,C>& sampleSet1,
																									 const Cluster<T,I,C>& sampleSet2,
																									 double stopAfter = MAXDOUBLE) const
			throw (ProximityException&);


		/**
		 * Check whether this measure is using representatives in
		 * calculating proximities between sets or in calculating the
		 * proximity between a sample and a set.
		 **/
		bool usesRepresentatives(void) const { return _bUseRepresentatives; }
		/**
		 * Set the use "representatives flag". The default value is true.
		 * See getProximity() for more info.
		 **/
		void setUseRepresentatives(bool use) { _bUseRepresentatives = use; }

		/**
		 * See whether this measure is symmetric.
		 **/
		bool isSymmetric(void) { return _bSymmetric; }

	protected:
		/**
		 * A flag that indicates whether this measure uses cluster
		 * representatives.
		 **/
		bool _bUseRepresentatives;
		/**
		 * A flag that indicates whether this measure is a symmetric
		 * measure, i.e. P(a,b) = P(b,a), where P(x,y) is the proximity
		 * between samples x and y.
		 **/
		bool _bSymmetric;
 	};

	template <class T>
	template <class I, class C> double ProximityMeasure<T>::getProximity(const Sample<T,I,C>& sample,
																																			 const Sample<T,I,C>& model,
																																			 double stopAfter) const
		throw (ProximityException&)
	{
		return getProximity(sample.featureVector(),model.featureVector(),stopAfter);
	}
	
	template <class T>
	template <class I, class C> double ProximityMeasure<T>::getProximity(const Sample<T,I,C>& sample,
																																			 const Cluster<T,I,C>& sampleSet,
																																			 double stopAfter) const
		throw (ProximityException&)
	{
		if (_bUseRepresentatives)
			return getProximity(sample,sampleSet.representative(),stopAfter);
		
		double minDist = stopAfter;
		for (int i=sampleSet.getLength();i--;)
			{
				double dist = getProximity(sample.featureVector(),sampleSet[i].featureVector(),minDist);
				if (dist < minDist)
					minDist = dist;
			}
		return minDist;
	}

	template <class T>
	template <class I, class C> double ProximityMeasure<T>::getProximity(const Cluster<T,I,C>& sampleSet1,
																																			 const Cluster<T,I,C>& sampleSet2,
																																			 double stopAfter) const
		throw (ProximityException&)
	{
		if (_bUseRepresentatives)
			return getProximity(sampleSet1.representative(),sampleSet2.representative(),stopAfter);

		double minDist = stopAfter;
		for (int i=sampleSet1.getLength();i--;)
			for (int j=sampleSet2.getLength();j--;)
				{
					double dist = getProximity(sampleSet1[i].featureVector(),sampleSet2[j].featureVector(),minDist);
					if (dist < minDist)
						minDist = dist;
				}
		return minDist;
	}

	/**
	 * A general modifier for proximity measures. This modifier may be
	 * used to perform calculations on proximity values in a general
	 * way. Stl functionals (plus,minus,multiplies,divides) are used to
	 * indicate the performed operations. ProximityModifiers can be used
	 * as any proximity measure - they just modify the value returned
	 * by their internal proximity measure in a certain way.<p>
	 *
	 * Example:<br>
	 * <pre>
	 * //Create a proximity modifier that adds 1.0 to the Euclidean distance between samples
	 * ProximityModifier&lt;plus&lt;double&gt;,double&gt; adder(new EuclideanDistance&lt;double&gt;, 1.0);
	 *
	 * //Create a proximity modifier that divides the Euclidean distance between samples by two
	 * ProximityModifier&lt;divides&lt;double&gt;,double&gt; adder(new EuclideanDistance&lt;double&gt;, 2.0);
	 * </pre>
	 **/
	template <class operation, class T> class ProximityModifier : public ProximityMeasure<T>
	{
	public:
		/**
		 * Create a new modifier that uses the given proximity measure.
		 * Note that if you use this constructor, the memory pointed by
		 * <i>measure</i> is automatically released upon the deletion of
		 * the modifier.
		 *
		 * @param measure a pointer to the actual measure to be used
		 * @param value the value to modify the measured proximities with
		 **/
		ProximityModifier(ProximityMeasure<T>* measure, double value) : _measure(measure), _dValue(value), _bDelete(true) {}
		/**
		 * Create a new modifier that uses the given proximity measure.
		 *
		 * @param measure a reference to the actual measure to be used
		 * @param value the value to modify the measured proximities with
		 **/
		ProximityModifier(ProximityMeasure<T>& measure, double value) : _measure(&measure), _dValue(value), _bDelete(false) {}

		virtual ~ProximityModifier() { if (_bDelete) delete _measure; }

		double getProximity(const util::List<T>& sample, const util::List<T>& model, double stopAfter) const throw (ProximityException&)
		{
			return _operation(_measure->getProximity(sample,model,stopAfter),_dValue);
		}
		
	private:
		ProximityMeasure<T>* _measure;
		double _dValue;
		bool _bDelete;
		operation _operation;
	};
	
	
	/**
	 * Adder is a wrap-up measure that just adds a constant value to the
	 * value provided by another measure. It may be used in scaling the
	 * measured values.<p>
	 * Example:<br>
	 * <pre>
	 * //Memory autorelease mode - no need to store a pointer to the
	 * //new EuclideanDiscance object
	 * ProximityAdder&lt;double&gt; adder(new EuclideanDistance&lt;double&gt;, -1.0);
	 *
	 * //No memory autorelease
	 * EuclideanDistance&lt;double&gt; dist;
	 * ProximityAdder&lt;double&gt; adder(dist, -1.0);
	 *
	 * //Absolutely wrong! The local variable would be automatically
	 * //deleted. Your program will most likely dump core.
	 * EuclideanDistance&lt;double&gt; dist;
	 * ProximityAdder&lt;double&gt; adder(&dist, -1.0); //One extra '&'.
	 * </pre>
	 * This measure comes in handy when you have multiple feature vectors
	 * in each sample and need to scale the proximities.
	 *
	 * @deprecated Use ProximityModifier instead
	 * @see ProximityModifier
	 * @see ProximityMultiplier
	 * @see ProximityMatrix
	 **/
	template <class T> class ProximityAdder : public ProximityMeasure<T>
	{
	public:
		/**
		 * Create a new Adder that uses the given proximity measure. Note
		 * that if you use this constructor, the memory pointed by
		 * <i>measure</i> is automatically released upon the deletion of
		 * the Adder.
		 * @param measure a pointer to the actual measure to be used
		 * @param value the value to be added to measured proximities
		 **/
		ProximityAdder(ProximityMeasure<T>* measure, double value) : _measure(measure), _dValue(value), _bDelete(true) {}
		/**
		 * Create a new Adder that uses the given proximity measure.
		 * @param measure a reference to the actual measure to be used
		 * @param value the value to be added to measured proximities
		 **/
		ProximityAdder(ProximityMeasure<T>& measure, double value) : _measure(&measure), _dValue(value), _bDelete(false) {}

		virtual ~ProximityAdder() { if (_bDelete) delete _measure; }

		double getProximity(const util::List<T>& sample, const util::List<T>& model, double stopAfter) const throw (ProximityException&)
		{
			return _measure->getProximity(sample,model,stopAfter) + _dValue;
		}
		
	private:
		ProximityMeasure<T>* _measure;
		double _dValue;
		bool _bDelete;
	};
	
	/**
	 * Multiplier is a wrap-up measure that multiplies the value
	 * provided by another measure by a constant value. It may be used
	 * in scaling the measured values.<p>
	 * Example:<br>
	 * <pre>
	 * MyProximityMeasure myProximity; //Inherited from ProximityMeasure&lt;float&gt;
	 *
	 * //Create a kNN classifier that uses MyProximityMeasure and scales the
	 * //proximity values to the range [0,1]. The values min and max can
	 * //be obtained for example from ProximityMatrix.
	 * kNNClassifier&lt;float&gt; knn(trainingSamples,
	 *                          new ProximityMultiplier&lt;float&gt;(new ProximityAdder&lt;float&gt;(myProximity,
	 *                                                                                  -min),
	 *                                                         1/(max-min)));
	 * </pre>
	 * @deprecated Use ProximityModifier instead
	 * @see ProximityModifier
	 * @see ProximityAdder
	 * @see ProximityMatrix
	 * @see Classifier
	 **/
	template <class T> class ProximityMultiplier : public ProximityMeasure<T>
	{
	public:
		/**
		 * Create a new Multiplier that uses the given proximity measure. Note
		 * that if you use this constructor, the memory pointed by
		 * <i>measure</i> is automatically released upon the deletion of
		 * the Multiplier.
		 * @param measure a pointer to the actual measure to be used
		 * @param value the value by which the measured proximities are to
		 *        be multiplied
		 **/
		ProximityMultiplier(ProximityMeasure<T>* measure, double value) : _measure(measure), _dValue(value), _bDelete(true) {}
		/**
		 * Create a new Multiplier that uses the given proximity measure.
		 * @param measure a pointer to the actual measure to be used
		 * @param value the value by which the measured proximities are to
		 *        be multiplied
		 **/
		ProximityMultiplier(ProximityMeasure<T>& measure, double value) : _measure(&measure), _dValue(value), _bDelete(false) {}

		virtual ~ProximityMultiplier() { if (_bDelete) delete _measure; }

		double getProximity(const util::List<T>& sample, const util::List<T>& model, double stopAfter) const throw (ProximityException&)
		{
			return _measure->getProximity(sample,model,stopAfter) * _dValue;
		}
		
	private:
		ProximityMeasure<T>* _measure;
		double _dValue;
		bool _bDelete;
	};

	/**
	 * The standard Euclidean distance as a proximity measure.
	 **/
	template <class T> class EuclideanDistance : public ProximityMeasure<T>
	{
	public:
		double getProximity(const util::List<T>& sample, const util::List<T>& model,
												double stopAfter = MAXDOUBLE) const throw (ProximityException&);
	};

	template <class T> double EuclideanDistance<T>::getProximity(const util::List<T>& lst1,
																															 const util::List<T>& lst2,
																															 double stopAfter) const
		throw (ProximityException&)
	{
		if (lst1.getLength() != lst2.getLength())
			throw ProximityException("EuclideanDistance::getProximity(const util::List<T>&, const util::List<T>&, double): Feature vectors differ in length.");

		double sum = 0.0, tmp, limit = (stopAfter == MAXDOUBLE) ? stopAfter : stopAfter*stopAfter;
		for (int i=0;i<lst1.getLength();i++)
			{
				tmp = (double)lst1[i] - (double)lst2[i];
				sum += tmp*tmp;
				if (sum > limit)
					return sum;
			}
		return sqrt(sum);
	}	

	/**
	 * Squared Euclidean distance as a proximity measure.
	 **/
	template <class T> class SquaredEuclidean : public ProximityMeasure<T>
	{
	public:
		double getProximity(const util::List<T>& sample, const util::List<T>& model,
												double stopAfter = MAXDOUBLE) const throw (ProximityException&);
	};

	template <class T> double SquaredEuclidean<T>::getProximity(const util::List<T>& lst1,
																															const util::List<T>& lst2,
																															double stopAfter) const
		throw (ProximityException&)
	{
		if (lst1.getLength() != lst2.getLength())
			throw ProximityException("SquaredEuclidean::getProximity(const util::List<T>&, const util::List<T>&, double): Feature vectors differ in length.");

		double sum = 0.0, tmp;
		for (int i=0;i<lst1.getLength();i++)
			{
				tmp = (double)lst1[i] - (double)lst2[i];
				sum += tmp*tmp;
				if (sum > stopAfter)
					return sum;
			}
		return sum;
	}	

	/**
	 * Cumlog is a log-likelihood proximity measure. It is defined as C
	 * = -sum<sub>i=1..N</sub>(S<sub>i</sub>ln(M<sub>i</sub>)), where S
	 * and M represent the sample and model distributions, respectively.
	 * N is the length of the distributions. If the normalization flag
	 * is set to true, M and S are normalized prior to distance
	 * calculation.<p>
	 *
	 * Cumlog is an asymmetric proximity measure, i.e. P(a,b) is
	 * generally different from P(b,a).
	 **/
	template <class T=double> class Cumlog : public ProximityMeasure<T>
	{
	public:
		/**
		 * Construct a new cumlog proximity measure. The measure is set to
		 * use cluster representatives.
		 *
		 * @param minValue the minimum value for a histogram bin (1.0
		 * means that zero bins are simply discarded because log(1.0)=0).
		 * @param normalized the normalization flag
		 **/
		Cumlog(double minValue=1e-8, bool normalized=true) :
			ProximityMeasure<T>(true,false),
			_dMinValue(minValue), _bNormalize(normalized), _matpLookup(NULL), _iMaxValue(0) {}

		/**
		 * Destroy the measure. If there is a look-up table, it will be
		 * destroyed as well.
		 **/
		~Cumlog() { delete _matpLookup; }
		
		/**
		 * Set the minimum value for a distribution entry that will be
		 * encountered. All zeros in model histograms will be replaced
		 * with this value in proximity calculations. (Because otherwise a
		 * logarithm of zero would be calculated.) The default value is
		 * 10<sup>-8</sup>. Setting this value is particularty important
		 * with sparse distributions. If your distribution is a histogram
		 * calculated from an matrix, you may want to set this value to
		 * 1/(rows*cols).
		 *
		 * @exception InvalidArgumentException& if <i>value</i> is smaller
		 *            than or equal to zero.
		 **/
		void setMinValue(double value) throw (InvalidArgumentException&)
		{
			if (value <= 0)
				throw InvalidArgumentException("Cumlog::setMinValue(double): value must be greater than zero.");
			_dMinValue = value;
		}
		
		/**
		 * Get the minimum value that is substituted to zero-valued
		 * feature vector entries.
		 **/
		double getMinValue(void) const { return _dMinValue; }
		/**
		 * Set the normalized state. If true, then all feature vectors
		 * given to this proximity measure are normalized prior to
		 * calculating the proximity. If false, then the feature vectors
		 * are assumed to be normalized a priori. In the latter case, the
		 * computational performance of the measure will be significantly
		 * higher. The default value is true.
		 **/
		void setNormalized(bool normalize) { _bNormalize = normalize; }
		/**
		 * Get the normalize state.
		 **/
		bool isNormalized(void) const { return _bNormalize; }

		double getProximity(const util::List<T>& sample, const util::List<T>& model,
												double stopAfter = MAXDOUBLE) const throw (ProximityException&);

		/**
		 * To speed up the calculation of the log-likelihood measure, one
		 * may think of using a look-up table instead of the on-line
		 * calculation. If you know that the values in your samples will
		 * never exceed a certain maximum value, then a look-up table of
		 * finite size can be created. For example, when creating a
		 * histogram out of 32x32 images, there will never be a value
		 * greater than 1024 (or smaller than 0). There are thus 1025
		 * distinct values, and the look-up table will contain 1025x1025
		 * entries. Alternatively, you may use the look-up table with
		 * normalized feature values. Each feature vector element is
		 * assumed to be in the range [0,1]. The values are multiplied by
		 * the given maximum value, rounded to the closest integer and
		 * used as a look-up table index.<p>
		 *
		 * Once created, the look-up table will be used for all subsequent
		 * proximity calculations. If you want to return to the on-line
		 * calculation mode, call releaseLookupTable(). Note that with
		 * normalized features, the proximity calculation is "quantized",
		 * resulting in somewhat more inaccurate proximity calculations.
		 *
		 * @param maxValue the maximum value a feature can attain.
		 **/
		void createLookupTable(int maxValue);

		/**
		 * Release the look-up table. All subsequent proximity
		 * calculations will be made on-line.
		 **/
		void releaseLookupTable() { delete _matpLookup; }

	private:
		double _dMinValue;
		bool _bNormalize;
		util::Matrix<double>* _matpLookup;
		int _iMaxValue;
	};

	template <class T> void Cumlog<T>::createLookupTable(int maxValue)
	{
		delete _matpLookup;

		int values = maxValue+1;
		_matpLookup = new util::Matrix<double>(values,values);
		double* ptr = _matpLookup->getData();
		double size = maxValue;
		for (int s=0;s<values;s++)
			for (int m=0;m<values;m++,ptr++)
				{
					double ds = (s > 0) ? s/size : _dMinValue;
					double dm = (m > 0) ? m/size : _dMinValue;
					*ptr = ds*log(dm);
				}
		_iMaxValue = maxValue;
	}

	template <class T> double Cumlog<T>::getProximity(const util::List<T>& sample,
																										const util::List<T>& model,
																										double stopAfter) const
		throw (ProximityException&)
	{
		if (sample.getLength() != model.getLength())
			throw ProximityException("Cumlog::getProximity(const util::List<T>&, const util::List<T>&, double): Feature vectors differ in length.");

		double sum = 0.0;

		if (!_bNormalize)
			{
				if (_matpLookup)
					{
						for (int i=0;i<sample.getLength();i++)
							{
								sum -= (*_matpLookup)(int(sample[i]*_iMaxValue + 0.5),int(model[i]*_iMaxValue + 0.5));
								if (sum > stopAfter)
									return sum;
							}
					}
				else
					{
						for (int i=0;i<sample.getLength();i++)
							{
								double s = double((sample[i] > _dMinValue) ? sample[i] : _dMinValue);
								double m = double((model[i] > _dMinValue) ? model[i] : _dMinValue);
								sum -= s*log(m);
								if (sum > stopAfter)
									return sum;
							}
					}
			}
		else
			{
				if (_matpLookup)
					{
						for (int i=0;i<sample.getLength();i++)
							{
								sum -= (*_matpLookup)(int(sample[i]),int(model[i]));
								if (sum > stopAfter)
									return sum;
							}
					}
				else
					{
						double modelSum = util::Math::sum(model), sampleSum = util::Math::sum(sample);
						double limit = (stopAfter == MAXDOUBLE) ? stopAfter : stopAfter * sampleSum;

						double tmp = _dMinValue * sampleSum;
						double tmp2 = _dMinValue * modelSum;
						for (int i=0;i<sample.getLength();i++)
							{
								double s = double((sample[i] > tmp) ? sample[i] : tmp);
								double m = double((model[i] > tmp2) ? model[i]/modelSum : _dMinValue);
								sum -= s*log(m);
								if (sum > limit)
									break;
							}
						sum /= sampleSum;
					}
			}
		return sum;
	}	
	

	/**
	 * JDDistance (Jeffrey's Divergence) is a statistical dissimilarity
	 * measure. It is defined as JD =
	 * -sum<sub>i=1..N</sub>(S<sub>i</sub>ln(2*S<sub>i</sub>/(M<sub>i</sub>+S<sub>i</sub>))
	 * +
	 * M<sub>i</sub>ln(2*M<sub>i</sub>/(M<sub>i</sub>+S<sub>i</sub>))),
	 * where S and M represent the sample and model distributions,
	 * respectively. N is the length of the distributions. If the
	 * normalization flag is set to true, M is normalized prior to
	 * distance calculation.
	 **/
	template <class T=double> class JDDistance : public ProximityMeasure<T>
	{
	public:
		/**
		 * Construct a new Jeffrey's Divergence proximity measure. The
		 * measure is normalized by default.
		 **/
		JDDistance() : _bNormalize(true), _dMinValue(1e-8) {}
		
		/**
		 * Set the minimum value for a distribution entry that will be
		 * encountered. All zeros in histograms will be replaced with this
		 * value in proximity calculations. (Because otherwise a logarithm
		 * of zero would be calculated.) Setting this value is
		 * particularty important with sparse distributions. If your
		 * distribution is a histogram calculated from an matrix, you may
		 * want to set this value to 1/(rows*cols).
		 *
		 * @exception InvalidArgumentException& if <i>value</i> is smaller
		 *            than or equal to zero.
		 **/
		void setMinValue(double value) throw (InvalidArgumentException&)
		{
			if (value <= 0)
				throw InvalidArgumentException("JDDistance::setMinValue(double): value must be greater than zero.");
			_dMinValue = value;
		}
		
		/**
		 * Get the minimum value that is substituted to zero-valued
		 * feature vector entries.
		 **/
		double getMinValue(void) const { return _dMinValue; }

		/**
		 * Set the normalized state. If true, then all feature vectors
		 * given to this proximity measure are normalized prior to
		 * calculating the proximity. If false, then the feature vectors
		 * are assumed to be normalized a priori. In the latter case, the
		 * computational performance of the measure will be significantly
		 * higher. The default value is true.
		 **/
		void setNormalized(bool normalize) { _bNormalize = normalize; }
		/**
		 * Get the normalize state.
		 **/
		bool isNormalized(void) { return _bNormalize; }
		double getProximity(const util::List<T>& sample, const util::List<T>& model,
												double stopAfter = MAXDOUBLE) const throw (ProximityException&);

	private:
		bool _bNormalize;
		double _dMinValue;
	};
	

	template <class T> double JDDistance<T>::getProximity(const util::List<T>& sample,
																												const util::List<T>& model,
																												double stopAfter) const
		throw (ProximityException&)
	{
		if (sample.getLength() != model.getLength())
			throw ProximityException("JDDistance::getProximity(const util::List<T>&, const util::List<T>&, double): Feature vectors differ in length.");

		double sum = 0.0;

		if (!_bNormalize)
			{
				for (int i=0;i<sample.getLength();i++)
					{
						double si = sample[i], mi = model[i];
						if (si == 0)
							si = _dMinValue;
						if (mi == 0)
							mi = _dMinValue;
						double denom = si+mi/2;
						
						sum += si*log(si/denom) + mi*log(mi/denom);
						if (sum > stopAfter)
							return sum;
					}
			}
		else
			{
				double modelSum = util::Math::sum(model), sampleSum = util::Math::sum(sample);

				for (int i=0;i<sample.getLength();i++)
					{
						double si = sample[i]/sampleSum, mi = model[i]/modelSum;
						if (si == 0)
							si = _dMinValue;
						if (mi == 0)
							mi = _dMinValue;
						double denom = si+mi/2;

						sum += si*log(si/denom) + mi*log(mi/denom);
						if (sum > stopAfter)
							return sum;
					}
			}
		return sum;
	}	
	
	/**
	 * Histogram intersection calculates the intersection between two
	 * feature distributions. The result of this measure is HI =
	 * 1-sum<sub>i=1..N</sub>(min(S<sub>i</sub>,M<sub>i</sub>)), where S
	 * and M represent the sample and model distributions, respectively.
	 * N is the length of the distributions.
	 **/
	template <class T=int> class HistogramIntersection : public ProximityMeasure<T>
	{
	public:
		/**
		 * Construct a new HistogramIntersection proximity measure. The
		 * measure is normalized by default.
		 **/
		HistogramIntersection() : _bNormalize(true) {}
		
		/**
		 * Set the normalized state. If true, then all feature vectors
		 * given to this proximity measure are normalized prior to
		 * calculating the proximity. If false, then the feature vectors
		 * are assumed to be normalized a priori. In the latter case, the
		 * computational performance of the measure will be significantly
		 * higher. The default value is true.
		 **/
		void setNormalized(bool normalize) { _bNormalize = normalize; }
		/**
		 * Get the normalize state.
		 **/
		bool isNormalized(void) { return _bNormalize; }
		double getProximity(const util::List<T>& sample, const util::List<T>& model,
												double stopAfter = MAXDOUBLE) const throw (ProximityException&);

	private:
		bool _bNormalize;
	};

	template <class T> double HistogramIntersection<T>::getProximity(const util::List<T>& sample,
																																	 const util::List<T>& model,
																																	 double stopAfter) const
		throw (ProximityException&)
	{
		if (sample.getLength() != model.getLength())
			throw ProximityException("HistogramIntersection::getProximity(const util::List<T>&, const util::List<T>&, double): Feature vectors differ in length.");

		double sum = 0.0;

		if (!_bNormalize)
			{
				for (int i=sample.getLength();i--;)
					sum += minimum(sample[i], model[i]);
			}
		else
			{
				double sampleSum = util::Math::sum(sample), modelSum = util::Math::sum(model);

				for (int i=sample.getLength();i--;)
					sum += minimum(sample[i]/sampleSum, model[i]/modelSum);
			}

		return 1-sum;
	}
}

#endif