/* -*- C++ -*-
*
* ---------------------------------------------------------------------
* $Id: teststatistics.cpp,v 1.2 2002/11/01 16:18:09 cag Exp $
* ---------------------------------------------------------------------
*
* Copyright (C) 2000-2002 Niv Drory <drory@usm.uni-muenchen.de>
*                         Claus A. Goessl <cag@usm.uni-muenchen.de>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2, or (at your option)
* any later version.
*
* This program 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, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA
*
* ---------------------------------------------------------------------
*
*/

#define LTL_RANGE_CHECKING

#include <cfloat>

#include <ltl/marray.h>
#include <ltl/statistics.h>

#include <iostream>


using namespace ltl;

using std::cout;
using std::endl;


int main(int argc, char **argv)
{
   cerr << "Testing MArray statistics (for now just average methods) ..." << endl;
   

   MArray<float,2> A(10,10);
   A = indexPos(A, 1) + 10 * (indexPos(A, 2) - 1);
   
   
   MArray<float,2> B(A.shape());
   B = A;
   B(Range(5,6), Range(5,6)) = 0;
   
   LTL_ASSERT_(  DBL_EPSILON >= fabs( (average(A) - 50.5) / 50.5 ),
                "average without nan failed" );

   LTL_ASSERT_( DBL_EPSILON >= fabs( (average(B, 0.0f) - 50.5) / 50.5 ),
                "average with nan present failed" );

   LTL_ASSERT_( DBL_EPSILON >= fabs( (average(B, -1.0f) - 48.48) / 48.48 ), 
                "average with nan but not present failed" );
   
   for(int k = 4; k > 1; --k)
   {
      double mean, sigma;
      double kappa = double(k);
      LTL_ASSERT_( 100 == kappa_sigma_average(A, kappa, &mean, &sigma), 
                   "kappa sigma without nan did not yield correct remaining elements" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (mean - 50.5) / 50.5), 
                   "kappa sigma without nan did not yield correct mean" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (sigma - 29.011491975882) / 29.011491975882), 
                   "kappa sigma without nan did not yield correct sigma" );
      LTL_ASSERT_( 96 == kappa_sigma_average(B, kappa, 0.0f, &mean, &sigma), 
                   "kappa sigma with nan present did not yield correct remaining elements" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (mean - 50.5) / 50.5), 
                   "kappa sigma with nan present did not yield correct mean" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (sigma - 29.598008467854) / 29.598008467854), 
                   "kappa sigma with nan present did not yield correct sigma" );
      LTL_ASSERT_( 100 == kappa_sigma_average(A, kappa, -1.0f, &mean, &sigma), 
                   "kappa sigma with nan but not present did not yield correct remaining elements" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (mean - 50.5) / 50.5), 
                   "kappa sigma with nan but not present did not yield correct mean" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (sigma - 29.011491975882) / 29.011491975882), 
                   "kappa sigma with nan but not present did not yield correct sigma" );
   }

   {
      double mean, sigma;
      LTL_ASSERT_( 2 == kappa_sigma_average(A, 1.0, &mean, &sigma), 
                   "kappa sigma without nan did not yield correct remaining elements" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (mean - 50.5) / 50.5), 
                   "kappa sigma without nan did not yield correct mean" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (sigma - M_SQRT1_2) / M_SQRT1_2), 
                   "kappa sigma without nan did not yield correct sigma" );
      LTL_ASSERT_( 2 == kappa_sigma_average(B, 1.0, 0.0f, &mean, &sigma), 
                   "kappa sigma with nan present did not yield correct remaining elements" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (mean - 50.5) / 50.5), 
                   "kappa sigma with nan present did not yield correct mean" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (sigma - M_SQRT1_2) / M_SQRT1_2), 
                   "kappa sigma with nan present did not yield correct sigma" );
      LTL_ASSERT_( 2 == kappa_sigma_average(A, 1.0, -1.0f, &mean, &sigma), 
                   "kappa sigma with nan but not present did not yield correct remaining elements" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (mean - 50.5) / 50.5), 
                   "kappa sigma with nan but not present did not yield correct mean" );
      LTL_ASSERT_( FLT_EPSILON >= fabs( (sigma - M_SQRT1_2) / M_SQRT1_2), 
                   "kappa sigma with nan but not present did not yield correct sigma" );
   }
}