TestMeanDoubleHelper.hh

#ifndef __TestMeanDoubleHelper_hh__
#define __TestMeanDoubleHelper_hh__

#include <Base/Math/Random.hh>
#include <Base/Math/Distribution.hh>

namespace BIAS {

  bool GenerateRandomData(Random &MyRandom,
                          const unsigned num_data, const unsigned num_outliers,
                          const double sigma3, const unsigned min_sample_size, 
                          std::vector<double> &data, 
                          std::vector<bool> &inliers, double &mean);
  
  bool CheckResult(const std::vector<bool> &est_inliers, const double &est_mean,
                   const std::vector<bool> &gt_inliers, const double &gt_mean,
                   const double &gt_sigma3);
  
  ///////////////////////////////////////////////////////////////////////
  // implementation
  //////////////////////////////////////////////////////////////////////


  bool GenerateRandomData(Random &MyRandom,
                          const unsigned num_data, const unsigned num_outliers,
                          const double sigma3, const unsigned min_sample_size,
                          std::vector<double> &data, 
                          std::vector<bool> &inliers, double &mean)
  {    
    BIASASSERT(num_outliers<num_data);
    if (num_data>=min_sample_size){
      // there is one check where not enough samples are generated
      BIASASSERT(num_outliers+min_sample_size<=num_data);
    }
    data.resize(num_data);
    inliers.resize(num_data);
    const double range = 20.;
    mean = MyRandom.GetUniformDistributed(-range, range);

    unsigned generated_inliers(0), generated_outliers(0);
    unsigned num_outliers_left = num_outliers;
    unsigned num_inliers_left = num_data - num_outliers;
    
    for (unsigned i=0; i<num_data; i++){
      BIASASSERT((num_data-i) == (num_outliers_left + num_inliers_left));

      unsigned inlier_outlier = 
        MyRandom.GetUniformDistributedInt(0u, 
                                          num_outliers_left + 
                                          num_inliers_left - 1u);

      bool is_inlier = (inlier_outlier<num_inliers_left);
      double val = 0.;

      if (is_inlier){
        // generate an inlier
        int count = 0, max_tries = 100;
        do {
          val = MyRandom.GetNormalDistributed(mean, sigma3/3.0);
        } while (fabs(val-mean)>=sigma3 && count<max_tries); 
        if (count==max_tries){
          val = mean + sigma3*0.5;
        }
        BIASASSERT(num_inliers_left>0);
        num_inliers_left--;
      } else {
        // generate an outlier
        int count = 0, max_tries = 100;
        do {
          val = MyRandom.GetUniformDistributed(-range, range);
          count++;
        } while (fabs(val-mean)<sigma3 && count<max_tries); 
        if (count==max_tries){
          val = mean * 4.0*sigma3;
        }
        BIASASSERT(num_outliers_left>0);
        num_outliers_left--;
      }
      data[i] = val;
      inliers[i] = is_inlier;
      //    cout << val<<"("<<((is_inlier)?("+"):("-"))<< ")  ";
      if (fabs(data[i]-mean)<sigma3){
        generated_inliers++;
        BIASASSERT(inliers[i]);
      } else {
        generated_outliers++;
        BIASASSERT(!inliers[i]);
      }
    }
    //cout << endl;
    BIASASSERT(num_inliers_left==0);
    BIASASSERT(num_outliers_left==0);
    BIASASSERT(generated_inliers+generated_outliers == num_data);
    BIASASSERT(generated_outliers == num_outliers);

    // check if the std. deviation of the outliers is smaller than the 
    // std. deviation of the inliers
    std::vector<double> inl, outl;
    for (unsigned i=0; i<data.size(); i++){
      if (inliers[i]){ 
        inl.push_back(data[i]);
      } else {
        outl.push_back(data[i]);
      }
    }
    if (!inl.empty() && !outl.empty()){
      Distribution dinl, doutl;
      dinl.AddData(inl);
      doutl.AddData(outl);
      static int num_tries = 0;
      const int max_num_tries = 10;
      if (dinl.GetSigma()*2.0 > doutl.GetSigma()){
        // try again
        if (num_tries<max_num_tries){
          num_tries++;
          return GenerateRandomData(MyRandom, num_data, num_outliers, sigma3,
                                    min_sample_size, data, inliers, mean);
        } else {
          std::cerr<<"error generating random data. This should never happen\n";
          return false;
        }
      }
    }
    return true;
  }


  bool CheckResult(const std::vector<bool> &est_inliers, const double &est_mean,
                   const std::vector<bool> &gt_inliers, const double &gt_mean,
                   const double &gt_sigma3)
  {
    if (fabs(est_mean - gt_mean)>gt_sigma3){
      std::cerr << "mean differs: \n   gt_mean: "<<gt_mean<<"\n  est_mean: "
                <<est_mean<<std::endl;
      return false;
    }
    if (est_inliers.size()!=gt_inliers.size()){
      std::cerr << "size of inliers mismatch: "<<std::endl;
      return false;
    }
    unsigned num_inlier_mismatches = 0;
    for (unsigned i=0; i<gt_inliers.size(); i++){
      if (gt_inliers[i]!=est_inliers[i]){
        num_inlier_mismatches++;
      }
    }
    double mismatch_ratio = 
      (double)num_inlier_mismatches / (double)gt_inliers.size();
    double max_inlier_mismatch_ratio = 0.1;
    if ((1./(double)gt_inliers.size())>max_inlier_mismatch_ratio){
      max_inlier_mismatch_ratio = 1./(double)gt_inliers.size();
    }
    if (mismatch_ratio>max_inlier_mismatch_ratio){
      std::cerr << "inlier mismtach ratio too big: "<<num_inlier_mismatches
                <<" / "<< gt_inliers.size()<< " = "<<mismatch_ratio <<std::endl;
      if (num_inlier_mismatches>5){
        return false;
      }
    }

    return true;
  }

} // namespace


#endif // __TestMeanDoubleHelper_hh__