ExamplePowell.cpp

example for solving a system of n non-linear equations with n unkowns using the powell hybrid method.

Author

MIP

/* 
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 ExamplePowell.cpp
    @relates Minpack
    @ingroup g_examples    
    @brief  example for solving a system of n non-linear equations with n unkowns using the
      powell hybrid method.
    @author MIP
*/


#include <Base/Common/BIASpragma.hh>

#include <MathAlgo/Minpack.hh>
#include <Base/Math/Random.hh>
#include <Base/Math/Vector.hh>
#include <Base/Math/Matrix3x3.hh>
#include <Base/Math/Vector3.hh>

using namespace BIAS;
using namespace std;

Vector<double> m1(3), m2(3);
int num;

int myfunc(void *p, int n, const double *x, double *fvec, int iflag)
{
  cerr << setw(3)<< num<<" : "<<x[0]<<"\t"<<x[1]<<"\t"<<x[2]<<" :: \t";
  fvec[0]=m2[0]-sin(x[0])*m1[0];
  fvec[1]=m2[1]-cos(x[1])*m1[1];
  fvec[2]=m2[2]-sin(x[2])*m1[2];
  cerr << fvec[0]<<"\t"<<fvec[1]<<"\t"<<fvec[3]<<endl;
  num++;
  return 0;
}



//int main(int argc, char *argv[])
int main()
{
  int res=0;
  Random ran;
  num=0;
  Matrix3x3<double> mat;
  Vector<double> groundtruth(3, "0.09 0.015 -0.02");
  Vector<double> InitialGuess(3, "0.15 0.02 -0.05");
  Vector<double> Result(3, "0.0 0.0 0.0");

  // create matches
  m1[0]=ran.GetUniformDistributed(-20.0, 20.0);
  m1[1]=ran.GetUniformDistributed(-20.0, 20.0);
  m1[2]=ran.GetUniformDistributed(-20.0, 20.0);
  //  mat.Fill(0.0);
  mat.SetZero();
  mat[0][0]=sin(groundtruth[0]);
  mat[1][1]=cos(groundtruth[1]);
  mat[2][2]=sin(groundtruth[2]);
  cout << "mat: "<<mat<<endl;
  m2=mat*(Vector3<double>)m1;
  cout <<m1<<"  <-->  "<<m2<<"\tmat * m1: "<<mat*(Vector3<double>)m1<<endl;
    
  if (Powell(myfunc,NULL, InitialGuess, Result)!=0){
    BIASERR("error in Powell ");
    return -1;
  } 

  cout <<  "ground truth        : "<<groundtruth
       <<"\nnon linear estimate : "<<Result<<endl;

  return res;
}