#ifndef STATWIDGET_HH_
#define STATWIDGET_HH_

/**
 * Project     : CASAVA
 * Module      : $RCSfile: StatWidget.hh,v $
 * @author     : Tony Cox
 * Copyright   : Copyright (c) Illumina 2008, 2009. All rights reserved.
 *
 ** This software is covered by the "Illumina Genome Analyzer Software
 ** License Agreement" and the "Illumina Source Code License Agreement",
 ** and certain third party copyright/licenses, and any user of this
 ** source file is bound by the terms therein (see accompanying files
 ** Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
 ** Illumina_Source_Code_License_Agreement.pdf and third party
 ** copyright/license notices).
 *
 */

#include <cmath>

#include <cstdio>

namespace ca {
namespace variance_detection {


/**
 * @class StatWidget
 *
 * @brief This is just a wrapper arround IndelFinder functionality.
 *
 * Accumulate mean and standard dev using a single pass formula
 * Uses the cancellation-friendly formulae on p.26 of
 * Higham, Accuracy & Stability of Numerical Algorithms
 * Variable names follow his
 * $Header $
 */
class StatWidget {
public:
	/** Constructor
	 */
	StatWidget() : M_(0),Q_(0),sum_(0),max_(0),min_(0),k_(0) {}

	virtual ~StatWidget();
	void update(double x) {
		k_++;
		if (k_==1) {
			max_=x;
			min_=x;

			M_=x;
		} else {
			if (x>max_)
				max_=x;
			if (x<min_)
				min_=x;

			// important to do M before Q as Q uses previous iterate of M
			Q_+=((k_-1)*pow((x-M_), 2))/k_;
			M_+=((x-M_)/k_);
		}
		sum_+=x;
	} // ~update

	void clear(void) {
		k_=0;
		sum_=0;
		Q_=0;
	} // ~clear

	int getNumSamples(void) const {
		return k_;
	} // ~getNumSamples(void) const

	double getMin(void) const {
		return min_;
	} // ~double getMin(void) const

	double getMax(void) const {
		return max_;
	} // ~double getMax(void) const

	double getMean(void) const {
		return ((k_ != 0) ? (sum_ / ((double)k_)) : 0.0);
	} // ~getMean(void) const

	double getSampleVariance(void) const { // s^2/n
		return ((k_ != 0) ? (Q_ / (double)k_) : 0.0);
	} // ~getSampleVariance(void) const

	double getPopulationVarianceEstimator(void) const { // s^2/(n-1)
		return ((k_ > 1) ? (Q_ / ((double)k_ - 1)) : 0.0);
	} // ~getPopulationVarianceEstimator(void) const

	double getStandardDeviation(void) const { // sqrt of variance
		return sqrt(getSampleVariance());
	} // ~getStandardDeviation(void) const

	double getSampleStandardDeviation(void) const { // sqrt of pve
		return sqrt(getPopulationVarianceEstimator());
	} // ~getSampleStandardDeviation(void) const

	double getStandardError(void) const {
		return ((k_ != 0)
                ? (getSampleStandardDeviation() / sqrt((double)k_))
                : 0.0);
	} // ~getStandardError(void) const

	void print(void) const {
		printf("entries=%d min=%f max=%f mean=%f sd=%f se=%f\n",
				getNumSamples(), getMin(), getMax(), getMean(),
				getStandardDeviation(), getStandardError());
	} // ~print(void) const

private:
	double M_;
	double Q_;
	double sum_;
	double max_;
	double min_;
	int k_;
};


}
} // end namespace casava{ namespace { applications
#endif /*STATWIDGET_HH_*/