MonteCarloTransform.cpp

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS 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 Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include <MathAlgo/MonteCarloTransform.hh>
#include <MathAlgo/SVD.hh>
#include <MathAlgo/Lapack.hh>
#include <Base/Common/CompareFloatingPoint.hh>
#include <Base/Debug/Exception.hh>
#include <Base/Math/Matrix.hh>

#include <fstream>

using namespace BIAS;
using namespace std;

MonteCarloTransform::
MonteCarloTransform()
  : UncertaintyTransformBase(), NumSamplesPerDimension_(100),
    UseQuasiRandomNumbers_(true), uniformDistribution_(false)
{
  NewDebugLevel("D_MCT_SAMPLES");
  NewDebugLevel("D_MCT_WRITE_SRC_SAMPLES");
  NewDebugLevel("D_MCT_WRITE_DST_SAMPLES");
  NewDebugLevel("D_MCT_NUM_SAMPLES");
  NewDebugLevel("D_MCT_EIGENVALVEC");
  //AddDebugLevel("D_MCT_SAMPLES");
}


MonteCarloTransform::
~MonteCarloTransform()
{}


int MonteCarloTransform::
Transform(const Vector<double>& src_mean, 
          const Matrix<double>& src_cov,
          std::vector<Vector<double> >& dst_samples) const
{
  int res = 0;
  vector<Vector<double> > src_samples;
  const int num = GetSamples_(src_mean, src_cov, src_samples);
  dst_samples.resize(num);
#ifdef BIAS_DEBUG
  ofstream of;
  if (DebugLevelIsSet("D_MCT_WRITE_DST_SAMPLES")){
    of.open("dst_samples.txt");
  }
#endif
#ifdef WRITE_DST_SAMPLES
  ofstream off("dst_samples_p4.txt");
#endif
  int num_errs = 0;
  for (int i=0; i<num; i++){
    if ((res = Transform_(src_samples[i], dst_samples[i-num_errs]))!=0){
      //BIASERR("error transforming "<<src_samples[i]<<"  error code: "<<res);
      num_errs++;
      continue;
    }
#ifdef BIAS_DEBUG
    if (DebugLevelIsSet("D_MCT_WRITE_DST_SAMPLES") && res==0){
      of << setw(5) << i << "\t";
      for (unsigned d=0; d<dst_samples[i-num_errs].Size(); d++){
        of << dst_samples[i-num_errs][d] <<"\t";
      }
      of << endl;
    }
#endif
#ifdef WRITE_DST_SAMPLES
    if (res==0){
      off << setw(5) << i << "\t";
      for (unsigned d=0; d<dst_samples[i-num_errs].Size(); d++){
        off << dst_samples[i-num_errs][d] <<"\t";
      }
      off << endl;
    }
#endif
  }

  if (num_errs!=0){
    BCDOUT(D_MCT_NUM_SAMPLES, num_errs<<" errors in transform functions, "
           "using less samples ("<<num-num_errs<<")\n");
    dst_samples.resize(num-num_errs);
  }
#ifdef BIAS_DEBUG
  if (DebugLevelIsSet("D_MCT_WRITE_DST_SAMPLES")){
    of.close();
  }
#endif    
#ifdef WRITE_DST_SAMPLES
  off.close();
#endif 
  return (num_errs==0)?(0):(-num_errs);
}


int MonteCarloTransform::
Transform(const Vector<double>& src_mean, const Matrix<double>& src_cov,
          Vector<double>& dst_mean, Matrix<double>& dst_cov) const
{
  if (src_cov.IsZero(numeric_limits<double>::epsilon())){
    Transform_(src_mean, dst_mean);
    const int dim = dst_mean.Size();
    dst_cov.newsize(dim, dim);
    dst_cov.SetZero();
  } else {
    vector<Vector<double> > dst_samples;
    int res = Transform(src_mean, src_cov, dst_samples);
    if (res!=0){
      BIASERR("error transforming samples. error code: "<<res);
      //return res;
    }
    const int num = (int)dst_samples.size();
    BIASASSERT(num>0);
    const int dim = (int)dst_samples[0].Size();
    dst_cov = Matrix<double>(dim, dim, MatrixZero);
    dst_mean.newsize(dim);
    dst_mean.SetZero();
    
    for (int i=0; i<num; i++){
      dst_mean += dst_samples[i];
    }
    BIASASSERT(num>0);
    dst_mean /= (double)num;
#ifdef BIAS_DEBUG
    // transform returned 0, so results MUST NOT contain any invalid numbers!
    for (int d=0; d<dim; d++) {
      INFNANCHECK(dst_mean[d]);
    }
#endif
    
    BIASASSERT(num>1);
    Vector<double> diff;
    for (int i=0; i<num; i++){
      diff = dst_samples[i] - dst_mean;
      dst_cov += diff.OuterProduct(diff);
    }
    dst_cov /= (double)(num-1);
#ifdef BIAS_DEBUG
    // now we check finally whether cov is valid
    for (int row=0; row<dst_cov.num_rows(); row++) {
      for (int col=0; col<dst_cov.num_cols(); col++) {
        INFNANCHECK(dst_cov[row][col]);
      }
    }
#endif
  }
  return 0;
}

//#define MC_NORMALIZE
int MonteCarloTransform::
GetSamples_(const Vector<double>& src_mean, 
            const Matrix<double>& src_cov,
            vector<Vector<double> >& src_samples) const
{
  static BIAS::Random randomizer;
  Vector<double> eigenval;
  vector<Vector<double> > eigenvec;
  int res;
#ifdef MC_NORMALIZE
  double factor = src_cov.Trace();
  Matrix<double> scaled_cov = src_cov / factor; 
  res = Upper_symmetric_eigenvalue_solve(scaled_cov, eigenval, eigenvec);
#else 
  res = Upper_symmetric_eigenvalue_solve(src_cov, eigenval, eigenvec);
#endif
  if (res!=0){
    BIASERR("error computing eigenvalues: "<<res);
    return -1;
  }
  BCDOUT(D_MCT_EIGENVALVEC, "eignevalues: "<<eigenval<<endl);
    
  const int dim = (int)src_mean.Size();
#ifdef BIAS_DEBUG
  //cout << "cov: "<<src_cov<<endl;
#endif
  for (int d=0; d<dim; d++){
    BCDOUT(D_MCT_EIGENVALVEC, d<<": "<<eigenvec[d]<<endl);
    // cope with numerical inaccuracies, covariance matrices must have
    // eignevalues >=0 
    if (eigenval[d]<0.0 && eigenval[d]>-1e-22){
      BCDOUT(D_MCT_EIGENVALVEC, "["<<d<<"]  numerical inaccurate eigenvalue: "
             <<eigenval[d]<<endl<< "corresp eigenvec.NormL2(): "
             <<eigenvec[d].NormL2()<<endl);
      eigenval[d] = 0.0;
    }
#ifdef BIAS_DEBUG
    if (eigenval[d]<0.0) { 
      cout << "["<<d<<"]  illegal eigenvalue: "<<eigenval[d]<<endl; 
      cout << "corresp eigenvec: "<<eigenvec[d]<<endl;
    }
#endif
    BIASASSERT(eigenval[d]>=0.0);
#ifdef MC_NORMALIZE
    eigenval[d] *= factor;
#endif
    eigenval[d] = sqrt(eigenval[d]);
#ifdef BIAS_DEBUG
    if (Equal(eigenvec[d].NormL2(), 0.0)) { 
      cout << "illegal eigenvector: "<<eigenvec[d]<<endl;}
#endif
    BIASASSERT(!Equal(eigenvec[d].NormL2(), 0.0));
    eigenvec[d] /= eigenvec[d].NormL2();
  }

#ifdef BIAS_DEBUG
  //cout << "src_cov: "<<src_cov<<endl;
  ofstream of;
  if (DebugLevelIsSet("D_MCT_WRITE_SRC_SAMPLES")){
    of.open("src_samples.txt");
  }
#endif
  double r;
  const int num =  NumSamplesPerDimension_ * dim;
  cout << "using: "<<num<<" samples\n";
  vector<Vector<double> > quasi_std_normal_distr;
  vector<double> normal_distr;


  if (UseQuasiRandomNumbers_){
    if (uniformDistribution_) { 
      BIASWARNONCE("uniform-distribution sampling is experimental"
                   ", multi-dim. case needs verification");
   
      BIAS::Vector<double> themin(dim, -1.0 * sqrt(3.0)),
        themax(dim, sqrt(3.0));
      res = randomizer.GetQuasiUniformDistributed(dim, themin, themax, num,
                                                   quasi_std_normal_distr);
      BIASASSERT(res==0);
    }  else {
      res = randomizer.GetQuasiNormalDistributed(dim, num, quasi_std_normal_distr);
      BIASASSERT(res==0);
    }
#ifdef BIAS_DEBUG
    // evaluate statistics of chosen samples:
    Vector<double> testmean;
    Matrix<double> testcov;
    randomizer.GetMeanAndCovariance(quasi_std_normal_distr,
                             testmean, testcov);
    BCDOUT(D_MCT_NUM_SAMPLES,"QuasiRandom with mean zero and cov=id computed "
           <<testmean<<" and "<<testcov);
#endif
  } else {
    if (uniformDistribution_){
      randomizer.GetUniformDistributed(-1.0 * sqrt(3.0), sqrt(3.0),
                                num*dim, normal_distr); 
    } else {
      randomizer.GetNormalDistributed(0.0, 1.0, num*dim, normal_distr);
    }
#ifdef BIAS_DEBUG
    // evaluate statistics of chosen samples:
    double testmean;
    double testcov;
    randomizer.GetMeanAndCovariance(normal_distr, testmean, testcov);
    BCDOUT(D_MCT_NUM_SAMPLES,"Random with mean zero and cov=id computed "
           <<testmean<<" and "<<testcov);
#endif
  }

  int jj=0;
  Vector<double> zero(dim, 0.0);
  src_samples.resize(num, zero);
  for (int i=0; i<num; i++){
    for (int d=0; d<dim; d++){

      if (UseQuasiRandomNumbers_){
        r = quasi_std_normal_distr[i][d] * eigenval[d];
      } else {
        //r = randomizer_.GetNormalDistributedGSL(0.0, eigenval[d]);
        r = normal_distr[jj++]*eigenval[d];
      }
/*
#ifdef BIAS_DEBUG
      // sanity checks
      if (src_samples[i].NormL2()>1e8){
        BIASERR("before adding vec: ["<<i<<"] "<<src_samples[i]<<" d = "<<d);
        BIASABORT; 
      }
#endif
*/
      src_samples[i] += r * eigenvec[d];
/*
#ifdef BIAS_DEBUG
      if (src_samples[i].NormL2()>1e8){
        BIASERR("before adding mean: random number "<<r<<"  "<<eigenval[d]
                <<"  "<<eigenvec[d]<<"  "<<src_samples[i]);
        BIASABORT; 
      }
#endif
*/
    }

    src_samples[i] += src_mean;
/*
#ifdef BIAS_DEBUG
    if (src_samples[i].NormL2()>1e8){
      BIASERR("after adding mean:   "<<src_samples[i]<<" "<<src_mean);
      BIASABORT; 
    }
#endif
*/

#ifdef BIAS_DEBUG
    if (DebugLevelIsSet("D_MCT_WRITE_SRC_SAMPLES")){
      for (int d=0; d<dim; d++){
        of << src_samples[i][d] <<"\t";
      }
      of << endl;
    }
#endif
    BCDOUT(D_MCT_NUM_SAMPLES, "x"<<flush);
  }
#ifdef BIAS_DEBUG
  if (DebugLevelIsSet("D_MCT_WRITE_SRC_SAMPLES")){
    of.close();
  }
#endif
  BCDOUT(D_MCT_NUM_SAMPLES, "\nusing "<<num<<" src_samples\n");
/*
#ifdef BIAS_DEBUG
  cout << "checking first two moments of distribution of "<<num<<" samples\n";
  // compute covariance and mean from samples, they should coincide with
  // src_mean and src_cov
  Vector<double> check_mean(dim, 0.0), diff;
  Matrix<double> check_cov(dim, dim, MatrixZero);
  for (int i=0; i<num; i++){
    //cout << i<<": "<<src_samples[i]<<endl;
    check_mean += src_samples[i];
  }
  check_mean /= (double)num;

  for (int i=0; i<num; i++){
    diff = src_samples[i] - check_mean;
    check_cov += diff.OuterProduct(diff);
  }
  check_cov /= (double)(num-1);

  diff = src_mean-check_mean;
  cout << "src_mean: "<<src_mean<<"\ncheck_mean: "<<check_mean
       << "diff: "<<diff<<"\t norm: "<<diff.NormL2()<<endl;
  cout << "src_cov: "<<src_cov
       << "check_cov: "<<check_cov<<endl
       << "src_cov - check_cov: "<<src_cov-check_cov<<endl;
#endif
*/
  return num;
}

void MonteCarloTransform::
Skew(const std::vector<Vector<double> >& dst_samples,
     double skewnness) const
{
  BIASERR("unfinished");
  BIASABORT;
}