ExampleUncertaintyTransform.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
*/

/** @example ExampleUncertaintyTransform.cpp
    @relates UnscentedTransform, MonteCarloTransform
    @ingroup g_examples    
    @brief  this is an example transforming uncertainty with UnscentedTransform
    @author woelk 07/2004 
*/

#include <MathAlgo/UnscentedTransform.hh>
#include <MathAlgo/MonteCarloTransform.hh>
#include <Base/Math/Random.hh>

using namespace BIAS;
using namespace std;

Matrix<double> T;

void MyConstTransform(const Vector<double>& src, Vector<double>& dst) 
{
  const int dim = src.Size();
  T.newsize(dim, dim);
  T.SetIdentity();
  dst = T * src;
}

void MyLinTransform(const Vector<double>& src, Vector<double>& dst) 
{
  const int dim = src.Size();
  T.newsize(dim, dim);
  const double c[] = { 1.0, 2.0, 0.0 };
  const int numc = 3;
  const int num = dim*dim;
  for (int i=0; i<num; i++){
    T.GetData()[i] = c[i%numc];
  }
  dst = T * src;
}

void MySqTransform(const Vector<double>& src, Vector<double>& dst) 
{
  const double c[] = { 1.57, 3.14, 6.28 };
  //const double c[] = { 1.0, 2.0, 0.0 };
  const int dim = src.Size();
  T.newsize(dim, dim);
  dst.newsize(dim);
  for (int i=0; i<dim; i++){
    dst[i] = c[0] + c[1] * src[i] + c[2] * src[i] * src[i];
  }
}

void MyTransform(const Vector<double>& src, Vector<double>& dst) 
{
  MyLinTransform(src, dst);
    //MyConstTransform(src, dst);
    //MySqTransform(src, dst);
}

/** \cond HIDDEN_SYMBOLS */
class MyUnscentedTransform
  : public UnscentedTransform
{
public:
  MyUnscentedTransform() : UnscentedTransform() { 
    /*AddDebugLevel(D_UT_TRANSFORM); 
    AddDebugLevel(D_UT_PARAM); 
    AddDebugLevel(D_UT_SIGMA_POINTS); */
  };

  virtual ~MyUnscentedTransform() {};

protected:

  virtual int Transform_(const Vector<double>& src,
                          Vector<double>& dst) const
  {
    MyTransform(src, dst);
    return 0;
  }
};

class MyMonteCarloTransform
  : public MonteCarloTransform
{
public:
  MyMonteCarloTransform() : MonteCarloTransform() { 
    /*AddDebugLevel(D_UT_TRANSFORM); 
    AddDebugLevel(D_UT_PARAM); 
    AddDebugLevel(D_UT_SIGMA_POINTS); */
  };

  virtual ~MyMonteCarloTransform() {};

protected:

  virtual int Transform_(const Vector<double>& src,
                          Vector<double>& dst) const
  {
    MyTransform(src, dst);
    return 0;
  }
};
/** \endcond */





int main(int /*argc*/, char ** /*argv*/)
{
  Vector<double> src, dst, gt_dst;
  Matrix<double> src_cov, dst_cov, gt_dst_cov;

  // generate src and src_cov using random number generator
  const int dim = 2;
  Random rand;
  src.newsize(dim);
  src_cov.newsize(dim, dim);
  for (int r=0; r<dim; r++){
    src[r] = rand.GetUniformDistributed(-10.0, 10.0);
    for (int c=r; c<dim; c++){
      src_cov[r][c] = src_cov[c][r] = 
        rand.GetUniformDistributed(-1.0, 1.0);
    }
    src_cov = src_cov * src_cov; 
  }
  /*src[0] = 1.0;
  src[1] = 2.0;
  src_cov[0][0] = 0.5;
  src_cov[1][1] = 0.75;
  src_cov[0][1] = src_cov[1][0] = 0.25; */
  

  cout << "src: "<<src<<endl;
  cout << "src_cov: "<<src_cov<<endl;
  
  MyTransform(src, gt_dst);
  gt_dst_cov = T * src_cov * T.Transpose();
  cout << "T: "<<T<<endl;
  cout << "gt: dst: "<<gt_dst<<endl;
  cout << "gt: dst_cov: "<<gt_dst_cov<<endl;

  MyUnscentedTransform mut;
  mut.SetBeta(2.0);
  mut.SetAlpha(0.001);

  if (mut.Transform(src, src_cov, dst, dst_cov)!=0){
    BIASERR("error transforming (unscented)");
    return -1;
  }
  cout << "ut: dst: "<<dst<<endl;
  cout << "ut: dst_cov: "<<dst_cov<<endl;

  cout << "diff: "<<(dst - gt_dst)<<endl;
  cout << "diff cov: "<<(dst_cov - gt_dst_cov)<<endl;

  MyMonteCarloTransform mmct;
  //mmct.AddDebugLevel("D_MCT_WRITE_SRC_SAMPLES");
  //mmct.AddDebugLevel("D_MCT_WRITE_DST_SAMPLES");
  //mmct.AddDebugLevel("D_MCT_EIGENVALVEC");
  mmct.SetNumSamplesPerDimension((int)1e2);
  mmct.SetUseQuasiRandomNumbers(true);
  if (mmct.Transform(src, src_cov, dst, dst_cov)!=0){
    BIASERR("error transforming (monte carlo)");
    return -1;
  }
  cout << "mc: dst: "<<dst<<endl;
  cout << "mc: dst_cov: "<<dst_cov<<endl;
  
  cout << "diff: "<<(dst - gt_dst)<<endl;
  cout << "diff cov: "<<(dst_cov - gt_dst_cov)<<endl;
  return 0;
}