MAth-Cloud-3.0.0/www.ipsocloud.com/intimamedia.com/cpp/MAth-Plaque-Dll-1.0.0/LibPlaque/Plaque/MeanEstimate.cpp

#include "MeanEstimate.h"

#include <cassert>
#include <cmath>
#include <limits>
#include <algorithm>


CDebugOutput::CDebugOutput() 
            : Singleton< CDebugOutput >(),
	      m_file( NULL ),
	      currentIndex( 0 ),
	      m_minIndex( -1 ),
	      m_maxIndex( -1 )
{
}

void CDebugOutput::Init()
{
  m_diameter.clear();
  m_estimate2.clear();
  m_distensibility.clear();
  m_distensibility20.clear();
  
  m_minIndex = -1;
  m_maxIndex = -1;
}

void CDebugOutput::OpenFile()
{
  m_file = fopen( "MathOutput.xls", "wt");
}


void CDebugOutput::WriteToFile()
{
  OpenFile();

  //ASSERT( m_distensibility20.GetSize() == m_diameter.GetSize() );
  assert( m_diameter.size() == m_estimate2.size() );
  
  int distSize = (int)m_distensibility20.size();

  double mean=0;
  {//mean
    int valide=0;
    double value=0;
    for(int i=0; i<distSize; i++)
    {
      std::map< int, double >::iterator it = m_distensibility.find( i );
      
      if(it != m_distensibility.end())
      {
        value= it->second;
        if(value > 0)
        {
          mean+=value;
          valide++;
        }
      }
    }
    mean /= valide;
  }

  double diam, esti, dist, dist20;

  for(int i=0; i<distSize; i++)
  {
          std::map< int, double >::iterator it = m_diameter.find( i );
  	  diam = esti = dist = dist20 = -10;
  	  if(it == m_diameter.end())
  	  {
  		  fprintf(m_file, "%d\t%f\t \t \t \t \tNo eim\n",i, mean);		  
  		  //fprintf(m_file,	  "%d\t%f\t \t \t \t \tNo eim\n",i, mean);		  
  	  }
  	  else
  	  {
	          diam = it->second;
		  std::map< int, double >::iterator it2 = m_estimate2.find( i );
		  esti = it2->second;
  		  if ( !CMeanEstimate::IsANumber( esti ))
  		  {
  			  fprintf(m_file, "%d\t%f\t%f\t \t \t \tNo eim within %lf\n",i, mean, diam, esti);
  		  } 
  		  else
  		  {	
		          std::map< int, double >::iterator itd = m_distensibility.find( i );
			  dist = itd->second;
			  itd = m_distensibility20.find( i );
			  dist20 = itd->second;
  			  if(dist20 == -2 )
  			  {
  				  fprintf(m_file, "%d\t%f\t%f\t%f\t%f\t \t OutOfBounds\n",i, mean, diam, esti, dist  );
  			  }
  			  else
  			  {
  				  fprintf(m_file, "%d\t%f\t%f\t%f\t%f\t%f\t%s\n",i, mean, diam, esti, dist, dist20, (m_minIndex==i? "MIN":(m_maxIndex==i? "MAX":"")) );
  			  }
  		  }
  	  }
  }

  fclose(m_file);
}



CMeanEstimate::CMeanEstimate()
{
}

CMeanEstimate::~CMeanEstimate(void)
{
}

void  CMeanEstimate::RemoveValues(std::vector< double > *A, double value)
{
  A->erase( remove_if( A->begin(), A->end(), bind2nd( std::equal_to< double >(), value ) ), A->end() );

/*
	double x;
	int size = (int)A->size();
	int count=0; //number removed;
	for(int i=0; i<size; i++)
	{
		x =  A->GetAt(i);
		if( x == value )
		{
			A->RemoveAt(i);
			count++;
			i--;
			size--;
		}
	}
*/
}

struct is_not_a_number : public std::unary_function< double, bool >
{
  bool operator() ( double x )
  {
    return !CMeanEstimate::IsANumber( x );
  }
};

void  CMeanEstimate::RemoveNanValues(std::vector< double > *A) // Nan = Not A Number
{
  A->erase( remove_if( A->begin(), A->end(), is_not_a_number() ), A->end() );

/*
	double x;
	int size = (int)A->GetSize();
	int count=0; //number removed;

	for(int i=0; i<size; i++)
	{
		x =  A->GetAt(i);
		if( ! CMeanEstimate::IsANumber(x) )
		{
			A->RemoveAt(i);
			count++;
			i--;
			size--;
		}
	}
*/
}


template < class Op1, class Op2, class Op3 >
class binary_compose : public std::unary_function< typename Op2::argument_type,
                                                   typename Op1::result_type >
{
  public:
  
    binary_compose( const Op1& o1, const Op2& o2, const Op3& o3 ) : f1( o1 ),
                                                                    f2( o2 ),
								    f3( o3 )
    {
    }
    
    typename Op1::result_type operator() ( const typename Op2::argument_type& x ) const
    {
      return f1( f2( x ), f3( x ) );
    }
    
  protected:
  
    Op1 f1;
    Op2 f2;
    Op3 f3;
};


template < class Op1, class Op2, class Op3 >
inline binary_compose< Op1, Op2, Op3 > compose2( const Op1& f1, const Op2& f2, const Op3& f3 )
{
  return binary_compose< Op1, Op2, Op3 >( f1, f2, f3 );
}


void  CMeanEstimate::RemoveOutOfBoundsValues(std::vector< double > *A, double percent)
{
	double mean = GetMeanEstimate(A, percent);
	double maxBound = mean * (1+percent);
	double minBound = mean * (1-percent);
	
  A->erase( remove_if( A->begin(), A->end(), 
                       compose2( std::logical_or< bool >(),
		                 std::bind2nd( std::less< double >(), minBound ),
		                 std::bind2nd( std::greater< double >(), maxBound ) ) ), 
		       A->end() );
	
/*
	double x=0;
	
	int size = (int)A->GetSize();
	int count=0; //number removed;
	
	for(int i=0; i<size; i++)
	{
		x =  A->GetAt(i);
		if( x > maxBound || x < minBound )
		{
			A->RemoveAt(i);
			count++;
			i--;
			size--;
		}
	}
*/
}

void  CMeanEstimate::FindMinMaxWithinBoundsValues(std::vector< double > *A, double percent, double &min, double &max)
{
	std::vector< double > A2(*A);
	RemoveValues(&A2, -1);
	RemoveValues(&A2, 0);
	RemoveNanValues(&A2);

//	min = std::numeric_limits<double>::max();
//	max = std::numeric_limits<double>::min();
	min = INT_MAX;
	max = INT_MIN;
	double mean = GetMeanEstimate(&A2, percent);
	double maxBound = mean * (1+percent);
	double minBound = mean * (1-percent);
	CDebugOutput::getInstance().m_maxIndex = -1;
	CDebugOutput::getInstance().m_minIndex = -1;
	double x=0;
	int i = 0;
	double value=-1;
	
	std::vector< double >::iterator 
	  it = A->begin(),
	  ie = A->end();
	
	while (it != ie)
	{
		x =  *it;
		std::map< int, double >::iterator iv = CDebugOutput::getInstance().m_distensibility.find( i );
		assert(iv != CDebugOutput::getInstance().m_distensibility.end());
		value = iv->second;
		assert( x == value );
		if(    !CMeanEstimate::IsANumber(x)
			|| x<=0){	

				CDebugOutput::getInstance().m_distensibility20.insert( std::make_pair(i, -1 ) );
		}else{	
			if (x > maxBound || x < minBound )
			{
				CDebugOutput::getInstance().m_distensibility20.insert( std::make_pair(i, -2) );
			}else{
				CDebugOutput::getInstance().m_distensibility20.insert( std::make_pair(i, x) );
				if( x > max ){ 
					CDebugOutput::getInstance().m_maxIndex = i;
					max = x;
				}else{			
					if(x < min ){
						CDebugOutput::getInstance().m_minIndex = i;
						min = x;
					}
				}
			}
		}
	  i++;
          ++it;
	}	
	assert(CDebugOutput::getInstance().m_maxIndex > 0 );
	assert(CDebugOutput::getInstance().m_minIndex > 0 );
}

void  CMeanEstimate::FindMinMax(std::vector< double > *A, double &min, double &max)
{
//	min = std::numeric_limits<double>::max();
//	max = std::numeric_limits<double>::min();
	min = INT_MAX;
	max = INT_MIN;
	//min = max = A->GetAt(0);
	double x=0;
	std::vector< double >::iterator
	  it = A->begin(),
	  ie = A->end();
	
	while ( it != ie )
	{
		x = *it;
		if( x > max ) max = x;
		if(x>0 && x < min ) min = x;
		++it;
	}
}

double CMeanEstimate::GetMeanEstimate(std::vector< double > *A)
{	
	return GetMeanEstimate(A, VARIATION_COEF);
}

void CMeanEstimate::startFileOutput(char *filename)
{
	CDebugOutput::getInstance().m_file = fopen(filename, "a");
	assert(CDebugOutput::getInstance().m_file);
}

void CMeanEstimate::titleFileOutput(char * /*title*/, std::vector< double > * /*A*/, std::vector< double > * /*B*/)
{
//	fprintf(m_file, "\n\n%s\n",title);
//	fprintf(m_file, "Sizes\n %d \t %d \n\n", A->GetSize(), B->GetSize());
//	fprintf(m_file, "Mean\n %f \t %f \n\n", CMeanEstimate::GetMean(A), CMeanEstimate::GetMean(B));
//	fprintf(m_file, "StandardDeviation\n %f \t %f \n\n", CMeanEstimate::GetStandardDeviation(A), CMeanEstimate::GetStandardDeviation(B));
}

void CMeanEstimate::titleFileOutput(std::vector< double > *A)
{
	fprintf(CDebugOutput::getInstance().m_file, "%f\t%f\t%d\t",GetMean(A),GetMeanEstimate(A), (int)A->size());
//	fprintf(m_file, "\n\n%s\n",title);
//	fprintf(m_file, "Sizes\n %d \t %d \n\n", A->GetSize(), B->GetSize());
//	fprintf(m_file, "Mean\n %f \t %f \n\n", CMeanEstimate::GetMean(A), CMeanEstimate::GetMean(B));
//	fprintf(m_file, "StandardDeviation\n %f \t %f \n\n", CMeanEstimate::GetStandardDeviation(A), CMeanEstimate::GetStandardDeviation(B));
}

void CMeanEstimate::writeFileOutput(char* text)
{
	fprintf(CDebugOutput::getInstance().m_file,text);
}

void CMeanEstimate::endFileOutput()
{
	fprintf(CDebugOutput::getInstance().m_file, "\n");
	fclose(CDebugOutput::getInstance().m_file);
}

double CMeanEstimate::GetMeanEstimate(std::vector< double > *A, double similarityCoef)
{
	double mean=GetMean(A);
	double x=0, sum = 0 ;
	int count=0;
	//double maxDeviation = mean * similarityCoef;
	std::vector< double >::iterator
	  it = A->begin(),
	  ie = A->end();
	
	while ( it != ie )
	{
		x = *it;
		//if(maxDeviation> fabs(x-mean)){
		if( x*similarityCoef > fabs(x-mean) )
		{
			sum+=x;
			count++;
		}
		else
		{
			int xxx=0; xxx++;
		}
		
		++it;
	}
//	VERIFY(count);
	return sum / (double)count;
}

double CMeanEstimate::GetMeanEstimateFileOutput(std::vector< double > *A, double similarityCoef)
{
	double mean=GetMean(A);
	double x=0, sum = 0 ;
	int count=0;
	//double maxDeviation = mean * similarityCoef;
	std::vector< double >::iterator
	  it = A->begin(),
	  ie = A->end();
	
	while ( it != ie )
	{
		x =  *it;
		//if(maxDeviation> fabs(x-mean)){
		if( x*similarityCoef > fabs(x-mean) )
		{
			sum+=x;
			count++;
		}
		
		++it;
	}
//#define OUTPUT 
#ifdef OUTPUT
	fprintf(m_file, "%d",count);
#endif
	//VERIFY(count);
	return sum / (double)count;
}

void CMeanEstimate::PrintMeanEstimate(std::vector< double > *A, std::vector< double > *B)
{
	PrintMeanEstimate(A,B, VARIATION_COEF);
}

void CMeanEstimate::PrintMeanEstimate(std::vector< double > *A, std::vector< double > *B, double similarityCoef)
{
	double mean =GetMean(A);
	double meanB=GetMean(B);
	double x=0;
	int count=0, countB=0;
	//fprintf(m_file, "MeanEstimate, coef:%f\n %f \t%f\n\n",similarityCoef, GetMeanEstimate(A,similarityCoef), GetMeanEstimate(B,similarityCoef));
	//fprintf(m_file, "AA\t exclu=0\t\tII\t exclu=0\n");
	//double maxDeviation = mean * similarityCoef;
	//double maxDeviationB = meanB * similarityCoef;
	
	assert( B->size() <= A->size() );
	std::vector< double >::iterator
	  ia = A->begin(),
	  ae = A->end(),
	  ib = B->begin();

	while ( ia != ae )
	{
		x =  *ia;
		//if(maxDeviation > fabs(x-mean) ){
		if ( x*similarityCoef > fabs(x-mean) )
		{
			//sum+=x;
			//fprintf(m_file, "%f\t1\t", x);
		}
		else
		{
			count++;
			//fprintf(m_file, "%f\t0\t", x);
		}
		x =  *ib;
		//fprintf(m_file, "\t");
		//if(maxDeviationB> fabs(x-meanB) ){
		if( x*similarityCoef > fabs(x-meanB) )
		{
			//sumB+=x;
			//fprintf(m_file, "%f\t1\t", x);
		}
		else
		{
			countB++;
			//fprintf(m_file, "%f\t0\t", x);
		}
		//fprintf(m_file, "\n");
		
		++ia;
		++ib;
	}
		//fprintf(m_file, "\nnombre de points exclus:\t%d\t\t\t%d\n",count, countB );
}

bool CMeanEstimate::IsANumber(double x)
{
//	return ((x < (double)(std::numeric_limits<double>::max())) && 
//		(x > (double)(std::numeric_limits<double>::min())));
	return ((x < (double)(INT_MAX)) && 
		(x > (double)(INT_MIN)));
}

double CMeanEstimate::GetMean(std::vector< double > *A)
{
	double sum=0;
	double x=0;
	std::vector< double >::iterator
	  it = A->begin(),
	  ie = A->end();

	while ( it != ie )
	{
		x = *it;
		if( !IsANumber(x)  )
		{
			int xxx=0; xxx++;
		} else
//		if( x == std::numeric_limits<double>::quiet_NaN() ){
//			int xxx=0; xxx++;
//		}
//		if( abs(x) == std::numeric_limits<double>::infinity() ){
//			int xxx=0; xxx++;
//		}
		sum+= x;
		
		++it;
	}
	return sum / (double)A->size();
}

double CMeanEstimate::GetStandardDeviation(std::vector< double > *A)
{
	return sqrt(GetVariance(A));
}

double CMeanEstimate::GetVariance(std::vector< double > *A)
{
	double mean=GetMean(A);
	double sum=0, diff=0;
	std::vector< double >::iterator
	  it = A->begin(),
	  ie = A->end();
	
	while ( it != ie )
	{
		diff = *it - mean;
		sum+= diff*diff;
		++it;
	}
	return sum / (double)A->size();
}

//void CMeanEstimate::OutputToFile(char* filename, CArray<double> *A, CArray<double> *B)
////implementation rapide et specifique est repectant peu le principde objet !!
//{
//	FILE* file = fopen(filename, "a");
//
//	if(file==0){
//		ASSERT(file);
//	}else{		
//		fprintf(file, "Mean\n %f \t %f \n", CMeanEstimate::GetMean(A), CMeanEstimate::GetMean(B));
//		fprintf(file, "StandardDeviation\n %f \t %f \n", CMeanEstimate::GetStandardDeviation(A), CMeanEstimate::GetStandardDeviation(B));
//		fprintf(file, "\nValues:\nAA \tII\n");
//		ASSERT(A->GetSize()==B->GetSize());
//		for(int i=A->GetSize()-1; i; i--){
//			fprintf(file, "%g \t %g \n", A->GetAt(i), B->GetAt(i));
//		}
//		fprintf(file, "\n\n");
//		fclose(file);
//	}
//}