#include <cmath>
#include <iostream>
#include <random>
#include <vector>

#include <cpExtensions/Algorithms/IterativeGaussianModelEstimator.h>

// -------------------------------------------------------------------------
const unsigned int Dim = 1;
typedef double TScalar;
typedef
cpExtensions::Algorithms::
IterativeGaussianModelEstimator< TScalar, Dim > TEstimator;

// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
  if( argc < 4 )
  {
    std::cerr << "Usage: " << argv[ 0 ] << " mean std samples" << std::endl;
    return( 1 );

  } // fi
  TScalar mean = std::atof( argv[ 1 ] );
  TScalar var = std::atof( argv[ 2 ] );
  unsigned int samples = std::atoi( argv[ 3 ] );
  var *= var;

  // Prepare estimator
  TEstimator::Pointer estimator = TEstimator::New( );
 
  // Generate numbers
  std::random_device r;
  std::seed_seq seed{ r( ), r( ), r( ), r( ), r( ), r( ), r( ), r( ) };
  std::mt19937 e( seed );
  std::normal_distribution< > dist( mean, std::sqrt( var ) );
  double local_mean = double( 0 );
  std::vector< double > data;
  for( unsigned int s = 0; s < samples; ++s )
  {
    double v = dist( e );
    estimator->AddSample( v );
    local_mean += v;
    data.push_back( v );

  } // rof
  local_mean /= double( samples );

  double local_var = double( 0 );
  for( auto d = data.begin( ); d != data.end( ); ++d )
    local_var += ( *d - local_mean ) * ( *d - local_mean );
  local_var /= double( samples - 1 );
  
  // Show results
  std::cout
    << "Mean: "
    << mean << " <-> "
    << estimator->GetMean( )[ 0 ] << " <-> "
    << local_mean
    << std::endl;
  std::cout
    << "Var: "
    << var << " <-> "
    << estimator->GetCovariance( )[ 0 ][ 0 ] << " <-> "
    << estimator->GetUnbiasedCovariance( )[ 0 ][ 0 ] << " <-> "
    << local_var
    << std::endl;

  std::cout << "--------------------------------------" << std::endl;
  for( unsigned int s = 0; s < 15; ++s )
  {
    double v = dist( e );
    double d = std::sqrt( estimator->SquaredMahalanobis( v ) );
    std::cout << "Distante to " << v << " is " << d << std::endl;

  } // rof

  return( 0 );
}

// eof - $RCSfile$