ExampleLaguerre.cpp

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

Copyright (C) 2003, 2004    (see file CONTACTS 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 ExampleLaguerre.cpp
    @relates Laguerre
    @ingroup g_examples    
    @brief   Example for usage of Laguerre Solver for Polynomials
    @author MIP
*/

#include "Base/Debug/Error.hh"
#include <MathAlgo/Laguerre.hh>
#include <Base/Common/CompareFloatingPoint.hh>
#include <iostream>
#include <algorithm>
#ifdef WIN32
#include <functional>
#endif


using namespace BIAS;
using namespace std;

double frand(double min = -10.0 , double max = 10.0)
{
  double tmin = min;
  double tmax = max;
  const unsigned int res = 10000000;
  unsigned int i=rand()%res;
  double scale  = ((double) i) / ((double) res);
  double diff   = tmax - tmin;
  return(tmin+scale*diff);
}

//bitmask looks like 00001111 
unsigned int GetNumberOfBinaryOnes(unsigned int val,
                                   unsigned int mask)
{
  unsigned int bitmask = mask;
  unsigned int temp    = val;
  unsigned int count   = 0;
  while((bitmask)!=0)  //dont use the full size of unsigned int
  {
    if((temp&1)==1)
    {
      count++;
    }
    temp    = temp>>1;
    bitmask = bitmask>>1;
  }
  return count;
}

void GetIndexListFromBinaryZeros(unsigned int val,
                                 unsigned int mask,
                                 std::vector<unsigned int>& indices)
{
  indices.clear();
  unsigned int bitmask = mask;
  unsigned int temp    = val;
  unsigned int count   = 0;
  while((bitmask)!=0)  //dont use the full size of unsigned int
  {
    if((temp&1)==0)
    {
      indices.push_back(count);
    }
    temp    = temp>>1;
    bitmask = bitmask>>1;
    count++;
  }
}

// create a polynomial from random roots
// idea: use bit pattern of numbers between 0 and (2^degree)-1
//       to decide which roots are multiplied and added to get the coefficients
void CreateRandomPolynomial(unsigned int degree,
                            vector< complex<double> >& coeffs,
                            vector< complex<double> >& roots)
{
  BIASASSERT(degree>0);
  roots.clear();
  coeffs.clear();
  std::vector<unsigned int> tempIDs;
  unsigned int coeffidx ;
  //create random roots ... (x-z0)(x-z1)(x-z2)...(x-zdegree)
  for(unsigned int r=0; r<degree; r++)
  {
    roots.push_back(frand());
  }
  coeffs.resize(degree+1,complex<double>(0,0));  //2
//   //create coeffs for random roots
  unsigned int maxiter = 1<<(degree); 
  unsigned int bitmask = maxiter - 1;   // 00111
  for(unsigned int c=0; c<maxiter; c++)
  {
    GetIndexListFromBinaryZeros(c,bitmask,tempIDs);
    coeffidx = GetNumberOfBinaryOnes(c,bitmask);
    complex<double> someroots(1,0);
    for(unsigned int r=0;r<tempIDs.size();r++)
    {
      someroots = someroots * (-roots[tempIDs[r]]);
    }
 //    cout<<c<<"  "<<"coeffidx: "<<coeffidx<<endl;
    BIASASSERT(coeffidx<coeffs.size());
    coeffs[coeffidx]= coeffs[coeffidx] + someroots;
  }

}

struct complex_comp : 
    public binary_function< complex<double>, complex<double>, bool> 
{
  bool operator()(complex<double> x, complex< double> y) 
     { 
       if(x.real()<=y.real())
       {
         if(x.real()==y.real())
         {
           return(x.imag()<y.imag());
         }
         return(true);
       }
       return(false);
     }
};


double MeanError(vector< complex<double> >& v1,
                 vector< complex<double> >& v2,
                 double& min, double& max)
{
  BIASASSERT(v1.size()==v2.size())
      
  sort(v1.begin(),v1.end(),complex_comp());
  sort(v2.begin(),v2.end(),complex_comp());
  
//   cout<<endl;
//   copy(v1.begin(),v1.end(), ostream_iterator<complex<double> >(cout, " "));
//   cout<<endl<<"--------------------------"<<endl;
//   copy(v2.begin(),v2.end(), ostream_iterator<complex<double> >(cout, " "));
//   
  double error =0 ;
  double tmax=0, tmin=10000000;
  for(unsigned int roo=0;roo<v1.size(); roo++)
  {
    double n = norm(v2[roo]-v1[roo]);
    if(n>tmax)tmax= n;
    if(n<tmin)tmin= n;
    error+= n; //sum deviation from 0.0
  }
  max=tmax;
  min=tmin;
  return (error/v1.size());
  
}

bool Equal(const complex<double>& a, 
           const complex<double>& b)
{
  return((Equal(a.real(),b.real(),10e-12))&&
        ((Equal(a.imag(),b.imag(),10e-12))));
}


int main(int argc, char *argv[])
{
  LaguerreSolver l;
  vector< complex<double> > realcoeff;
  vector< complex<double> > expectedroots;
  vector< complex<double> > roots;
  
  unsigned int deg = 8;  //<------------------------------- hier!!!
  
  double max = -10e6; // small number
  double min = 10e6; // large number
  double meanmeanerror=0;
  double meanmin = 10000000;
  double meanmax = 0;
  double tmin;
  double tmax;
  for(unsigned run=0; run<100; run++)
  {
    CreateRandomPolynomial(deg,realcoeff,expectedroots);
    l.Solve(realcoeff, roots);
    double err;
    if(run==0)
      err=MeanError(expectedroots,roots,min,max);
    else
    {
      err=MeanError(expectedroots,roots,tmin,tmax);
     if(min>tmin)min=tmin;
     if(max<tmax)max=tmax;
    }
    cout<<err<<" min: "<< min<<" max: "<< max<<endl;
   
    if(err<meanmin)meanmin=err;
    if(err>meanmax)meanmax=err;
    meanmeanerror +=err;
  }
  meanmeanerror /= 100.0;
  cout<<"-----------------------------------------------------------"<<endl;
  cout<<"min error (no mean):"<< min <<endl;
  cout<<"max error (no mean):"<< max <<endl;
  cout<<"mean min error:"<< meanmin <<endl;
  cout<<"mean max error:"<< meanmax <<endl;
  cout<<"mean mean error:"<< meanmeanerror <<endl;
  return 0;
}