/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2014 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once
#if !defined(__MITSUBA_CORE_CHISQUARE_H_)
#define __MITSUBA_CORE_CHISQUARE_H_

#include <mitsuba/render/common.h>
#include <boost/tuple/tuple.hpp>
#include <boost/function.hpp>

MTS_NAMESPACE_BEGIN

/// Minimum expected cell frequency. Cells below this value will be pooled
#define CHISQR_MIN_EXP_FREQUENCY 5

/**
 * \brief Chi-square goodness-of-fit test on the sphere
 *
 * This class performs a chi-square goodness-of-fit test of the null hypothesis
 * that a specified sampling procedure produces samples that are distributed
 * according to a supplied density function. This is very useful to verify BRDF
 * and phase function sampling codes for their correctness. Currently, it
 * supports both 2D and discrete sampling methods and mixtures thereof.
 *
 * This implementation works by generating a large batch of samples, which are
 * then accumulated into rectangular bins in spherical coordinates. To obtain
 * reference bin counts, the provided density function is numerically
 * integrated over the area of each bin. Comparing the actual and reference
 * bin counts yields the desired test statistic.
 *
 * Given a probability distribution with the following interface
 *
 * \code
 * class MyDistribution {
 *     // Sample a (optionally weighted) direction. A non-unity weight
 *     // in the return value is needed when the sampling distribution
 *     // doesn't exactly match the implementation in pdf()
 *     boost::tuple<Vector, Float, EMeasure> generateSample() const;
 *
 *     /// Compute the probability density for the specified direction and measure
 *     Float pdf(const Vector &direction, EMeasure) const;
 * };
 * \endcode
 *
 * the code in this class might be used as follows
 *
 * \code
 * MyDistribution myDistrInstance;
 * ChiSquare chiSqr;
 *
 * // Initialize the tables used by the chi-square test
 * chiSqr.fill(
 *    boost::bind(&MyDistribution::generateSample, myDistrInstance),
 *    boost::bind(&MyDistribution::pdf, myDistrInstance, _1, _2)
 * );
 *
 * // Optional: dump the tables to a MATLAB file for external analysis
 * chiSqr.dumpTables("debug.m");
 *
 * if (!chiSqr.runTest())
 *    Log(EError, "Uh oh -- test failed, the implementation is probably incorrect!");
 * \endcode
 * \ingroup libcore
 */
class MTS_EXPORT_CORE ChiSquare : public Object {
public:
	/// Possible outcomes in \ref runTest()
	enum ETestResult {
		/// The null hypothesis was rejected
		EReject = 0,
		/// The null hypothesis was accepted
		EAccept = 1,
		/// The degrees of freedom were too low
		ELowDoF = 2
	};

	/**
	 * \brief Create a new Chi-square test instance with the given
	 * resolution and sample count
	 *
	 * \param thetaBins
	 *    Number of bins wrt. latitude. The default is 10
	 *
	 * \param phiBins
	 *    Number of bins wrt. azimuth. The default is to use
	 *    twice the number of \c thetaBins
	 *
	 * \param numTests
	 *    Number of independent tests that will be performed. This
	 *    is used to compute the Sidak-correction factor.
	 *
	 * \param sampleCount
	 *    Number of samples to be used when computing the bin
	 *    values. The default is \c thetaBins*phiBins*5000
	 */
	ChiSquare(int thetaBins = 10, int phiBins = 0,
			int numTests = 1, size_t sampleCount = 0);

	/// Get the log level
	inline ELogLevel getLogLevel() const { return m_logLevel; }

	/// Set the log level
	inline void setLogLevel(ELogLevel logLevel) { m_logLevel = logLevel; }

	/**
	 * \brief Set the tolerance threshold for bins with very low
	 * aggregate probabilities
	 *
	 * When the Chi-square test integrates the supplied probability
	 * density function over the support of a bin and determines that
	 * the aggregate bin probability is zero, the test would ordinarily
	 * fail if as much as one sample is placed in that bin in the
	 * subsequent sampling step. However, due to various numerical
	 * errors in a system based on finite-precision arithmetic, it
	 * may be a good idea to tolerate at least a few samples without
	 * immediately rejecting the null hypothesis. This parameter
	 * sets this threshold. The default value is \c number-of-samples*1e-4f
	 */
	inline void setTolerance(Float tolerance) { m_tolerance = tolerance; }

	/**
	 * \brief Fill the actual and reference bin counts
	 *
	 * Please see the class documentation for a description
	 * on how to invoke this function
	 */
	void fill(
		const boost::function<boost::tuple<Vector, Float, EMeasure>()> &sampleFn,
		const boost::function<Float (const Vector &, EMeasure)> &pdfFn);

	/**
	 * \brief Dump the bin counts to a file using MATLAB format
	 */
	void dumpTables(const fs::path &filename);

	/**
	 * \brief Perform the actual chi-square test
	 *
	 * \param pvalThresh
	 *     The implementation will reject the null hypothesis
	 *     when the computed p-value lies below this parameter
	 *     (default: 0.01f)
	 *
	 * \return A status value of type \ref ETestResult
	 */
	ETestResult runTest(Float pvalThresh = 0.01f);

	MTS_DECLARE_CLASS()
protected:
	/// Release all memory
	virtual ~ChiSquare();

	/// Functor to evaluate the pdf values in parallel using OpenMP
	static void integrand(
		const boost::function<Float (const Vector &, EMeasure)> &pdfFn,
			size_t nPts, const Float *in, Float *out) {
		#if defined(MTS_OPENMP)
		#pragma omp parallel for
		#endif
		for (int i=0; i<(int) nPts; ++i)
			out[i] = pdfFn(sphericalDirection(in[2*i], in[2*i+1]), ESolidAngle)
				* std::sin(in[2*i]);
	}
private:
	ELogLevel m_logLevel;
	Float m_tolerance;
	int m_thetaBins, m_phiBins;
	int m_numTests;
	size_t m_sampleCount;
	Float *m_table;
	Float *m_refTable;
};

MTS_NAMESPACE_END

#endif /* __MITSUBA_CORE_CHISQUARE_H_ */
