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

#ifndef _GAUSSIAN_H
#define _GAUSSIAN_H

#include <Matrix.h>
#include <math.h>

namespace prapi
{
	/**
	 * A class that contains methods for calculating values from
	 * Gaussian distributions and for creating various forms of a
	 * Gaussian function.
	 **/
	class Gaussian
	{
	public:
		/**
		 * Create a simple, center-symmetric 2D Gaussian function. The
		 * center of the Gaussian "hill" is at the center of the returned
		 * matrix. The width of the Gaussian function can be controlled by
		 * two parameters: radius and width. Radius specifies the range of
		 * values the Gaussian is calculated on, and width specifies the
		 * width of the Gaussian itself. The default width is one, which
		 * results in a Gaussian function with a variance of 0.5. The size
		 * parameter effectively specifies the sampling resolution of the
		 * Gaussian function. The maximum value of the returned matrix is
		 * 1.
		 *
		 * @param radius calculate values between [-radius,radius]
		 * @param size the size of the output matrix (specifies sampling frequency)
		 * @param width the width of the gaussian (2*std<sup>2</sup>)
		 * @return a size x size matrix
		 **/
		static util::Matrix<double> create2DGaussian(double radius, int size, double width=1);

		/**
		 * Create a 2D normal distribution. The larger radius and matrix
		 * size you use, the closer the sum of the result will be to
		 * unity.
		 *
		 * @param radius calculate values between [-radius,radius]
		 * @param size of the output matrix
		 * @param deviation the standard deviation of the distribution
		 **/
		static util::Matrix<double> create2DNormal(double radius, int size, double deviation);
		
		/**
		 * Create a list of values containing a sampled 1-D representation
		 * of a Gaussian function. The first item of the returned list
		 * equals to G(-range), and the last to G(range). The size
		 * parameter specified the number of sample taken between (and
		 * including) these extrema.
		 *
		 * @see #create2DGaussian(double,int,double)
		 **/
		static util::List<double> create1DGaussian(double range, int size, double width=1);

		/**
		 * Get the value of the normal distribution N(mean,deviation) at
		 * x. If mean is set to 0 and deviation to 1, the distribution is
		 * the "standard" normal distribution N(0,1).
		 **/
		static double gaussian(double x, double mean=0, double deviation=1)
		{
			return M_2_SQRTPI / (2*M_SQRT2) / deviation * exp(-(x-mean)*(x-mean)/(2*deviation*deviation));
		}
		
		/**
		 * The "error function" (erf). Erf is encountered in integrating
		 * the normal distribution (normalized Gaussian). erf(x) = int
		 * -Inf->x G(t)dt, where G(t) is a Gaussian distribution with zero
		 * mean and unit variance. The value of the error function is
		 * calculated as a Maclaurin series, which may take quite a lot of
		 * time.
		 **/
		static double erf(double x);
	};
}

#endif