//------------------------------------------------------------------------------
//	Copyright (C) 2009 Humboldt-Universitaet zu Berlin
//
//	This library is free software; you can redistribute it and/or
//	modify it under the terms of the GNU Lesser General Public
//	License as published by the Free Software Foundation; either
//	version 2.1 of the License, or (at your option) any later version.
//
//	This library 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
//	Lesser General Public License for more details.
//
//	You should have received a copy of the GNU Lesser General Public
//	License along with this library; if not, write to the Free Software
//	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//------------------------------------------------------------------------------

/**
 * @file Regression.cpp
 * @author Ronald Kluth
 * @date created at 2009/03/01
 * @brief Implementation of odemx::statistics::Regression
 * @sa Regression.h
 * @since 3.0
 */

#include <odemx/statistics/Regression.h>
#include <odemx/data/ReportTable.h>
#include <odemx/data/Report.h>

#include <cmath>

namespace odemx {
namespace statistics {

//------------------------------------------------------construction/destruction

Regression::Regression( base::Simulation& sim, const data::Label& label )
:	Tab( sim, label )
,	sumX_( 0.0 )
,	sumY_( 0.0 )
,	sumSqX_( 0.0 )
,	sumXMulY_( 0.0 )
,	sumSqY_( 0.0 )
{
}

Regression::~Regression()
{
}

//---------------------------------------------------------------------modifiers

void Regression::update( double valueX, double valueY)
{
	Tab::update();

	sumX_ += valueX;
	sumY_ += valueY;
	sumSqX_ += std::pow( valueX, 2 );
	sumSqY_ += std::pow( valueY, 2 );
	sumXMulY_ += valueX * valueY;
}

void Regression::reset( base::SimTime time )
{
	Tab::reset( time );
	sumX_ = sumY_ = sumSqX_ = sumXMulY_ = sumSqY_ = 0.0;
}

//-----------------------------------------------------------------------getters

bool Regression::hasSufficientData() const
{
	return getUpdateCount() > 5;
}

double Regression::getXMean() const
{
	return getUpdateCount() ? sumX_ / getUpdateCount() : 0.0;
}

double Regression::getYMean() const
{
	return getUpdateCount() ? sumY_ / getUpdateCount() : 0.0;
}

double Regression::getDx() const
{
	return std::fabs( getUpdateCount() * sumSqX_ - std::pow( sumX_, 2 ) );
}

double Regression::getDy() const
{
	return std::fabs( getUpdateCount() * sumSqY_ - std::pow( sumY_, 2 ) );
}

double Regression::getResidualStandardDeviation() const
{
	// RES.ST.DEV
	double sd = std::sqrt( ( sumSqY_ - getIntercept() * sumY_
			- getEstimatedRegressionCoefficient() * sumXMulY_ )
			/ (getUpdateCount() - 2 ) );
	return sd;
}

double Regression::getEstimatedRegressionCoefficient() const
{
	// EST.REG.COEFF
	double a1 = ( getUpdateCount() * sumXMulY_ - sumX_ * sumY_ ) / getDx();
	return a1;
}

double Regression::getIntercept() const
{
	// INTERCEPT
	double a0 = ( sumY_ * sumSqX_ - sumX_ * sumXMulY_ ) / getDx();
	return a0;
}

double Regression::getRegressionCoefficientStandardDeviation() const
{
	// ST.DEV.REG.COEFF
	return ( getUpdateCount() * getResidualStandardDeviation()
			/ std::sqrt( ( getUpdateCount() - 2 ) * getDx() ) );
}

double Regression::getCorrelationCoefficient() const
{
	// CORR.COEFF
	double r2 = ( ( getUpdateCount() * sumXMulY_ - sumX_ * sumY_ )
			* ( getUpdateCount() * sumXMulY_ - sumX_ * sumY_ ) )
			/ ( getDx() * getDy() );
	return r2;
}

//------------------------------------------------------------------------report

void Regression::report( data::Report& report )
{
	using namespace data;
	ReportTable::Definition def;
	def.addColumn( "Name", ReportTable::Column::STRING );
	def.addColumn( "Reset at", ReportTable::Column::SIMTIME );
	def.addColumn( "Uses", ReportTable::Column::INTEGER );
	def.addColumn( "XBar", ReportTable::Column::REAL );
	def.addColumn( "YBar", ReportTable::Column::REAL );
	def.addColumn( "Dx", ReportTable::Column::REAL );
	def.addColumn( "Dy", ReportTable::Column::REAL );
	def.addColumn( "Residual StDev", ReportTable::Column::REAL );
	def.addColumn( "Est Reg Coeff", ReportTable::Column::REAL );
	def.addColumn( "Intercept", ReportTable::Column::REAL );
	def.addColumn( "Reg Coeff StDev", ReportTable::Column::REAL );
	def.addColumn( "Correlation Coeff", ReportTable::Column::REAL );
	ReportTable& table = report.getTable( "Regression Statistics", def );
	table
	<< getLabel()
	<< getResetTime()
	<< getUpdateCount()
	<< getXMean()
	<< getYMean()
	<< getDx()
	<< getDy()
	<< getResidualStandardDeviation()
	<< getEstimatedRegressionCoefficient()
	<< getIntercept()
	<< getRegressionCoefficientStandardDeviation()
	<< getCorrelationCoefficient();
}

} } // namespace odemx::statistics
