/** 
* @file MAPMultinomialDirichlet.h
// refactor-nice.pl: check this substitution
// old: * @brief map estimation of a multinomial using a dirichlet prior
* @brief std::map estimation of a multinomial using a dirichlet prior
* @author Erik Rodner
* @date 10/30/2008

*/
#include <assert.h>

#include "MAPMultinomialDirichlet.h"

using namespace OBJREC;
using namespace NICE;
using namespace std;

void MAPMultinomialDirichlet::estimate ( NICE::Vector & mapEstimate,
		    const VVector & likelihoodDistributionSamples,
		    const VVector & priorDistributionSamples,
		    double priorInfluence )
{
    assert ( likelihoodDistributionSamples.size() == 1 );

    priorInfluence = 1.0 / priorInfluence;

    const NICE::Vector & mlEstimate = likelihoodDistributionSamples[0];

    double mlEstimateSum = 0.0;
    for ( uint k = 0 ; k < (uint)mlEstimate.size(); k++ )
	mlEstimateSum += mlEstimate[k];

    NICE::Vector mu;
    for ( VVector::const_iterator i = priorDistributionSamples.begin();
			    i != priorDistributionSamples.end(); i++ )
    {
	const NICE::Vector & x = *i;
	if ( mu.size() == 0 ) mu = x;
	else                  mu = mu + x;
    }
    mu = mu * (priorInfluence/priorDistributionSamples.size());

    // ------ CODE SAFETY
    double muSum = 0.0;
    for ( uint k = 0 ; k < (uint)mu.size(); k++ )
	muSum += mu[k];
    assert ( fabs(muSum - priorInfluence) < 10e-9 );
    // ------ END CODE SAFETY

    // mu is a rough estimate of the parameter vector alpha of a dirichlet distribution

    mapEstimate.resize(mlEstimate.size());
    mapEstimate.set(0);

    double scale = mlEstimateSum + priorInfluence - mlEstimate.size();
    assert ( fabs(scale) > 10e-8 );

    for ( uint k = 0 ; k < (uint)mapEstimate.size() ; k++ )
	mapEstimate[k] = (mlEstimate[k] + mu[k] - 1) / scale;

}