ExamplePolynomialSolve.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 ExamplePolynomialSolve.cpp
    @relates PolynomialSolve
    @ingroup g_examples    
    @brief  example for PolynomialSolve usage
    @author MIP
*/


#include <Base/Common/BIASpragma.hh>

#include <MathAlgo/PolynomialSolve.hh>
#include <Base/Math/Random.hh>
#include <algorithm>

using namespace BIAS;
using namespace std;

bool debug=false;
double maxerr=1e-8;

namespace std{
  bool operator<(complex<double>& left, complex<double>& right)
  {
    return (left.real()<right.real());
  }
  bool operator<(const complex<double>& left, const complex<double>& right)
  {
    return (left.real()<right.real());
  }
}

void solout(vector<POLYNOMIALSOLVE_TYPE>& sol)
{
  sort(sol.begin(), sol.end());
  int s=(int)sol.size();
  cout.precision(20);
  for (int i=0; i<s; i++)
    if (debug) cout << sol[i] << "  ";
  if (debug) cout << endl;
}

void solout(vector<complex<POLYNOMIALSOLVE_TYPE> >& sol)
{
  sort(sol.begin(), sol.end());
  int s=(int)sol.size();
  cout.precision(20);
  for (int i=0; i<s; i++)
    if (debug) cout << sol[i] << "  ";
  if (debug) cout << endl;
}


void errout(vector<POLYNOMIALSOLVE_TYPE>& sol, 
            vector<POLYNOMIALSOLVE_TYPE>& gtsol)
{
  sort(sol.begin(), sol.end());
  int s=(int)sol.size();
  cout.precision(20);
  if (debug) cout << "error: ";
  for (int i=0; i<s; i++){
    
   if (debug) cout << fabs(sol[i]-gtsol[i]) << "  ";
    if (fabs(sol[i]-gtsol[i])>maxerr) {
      cerr << "\n\n\nerror bigger than "<<maxerr<<"\n" <<endl; 
      debug=true;
      cout << "ground truth: ";
      solout(gtsol);
      cout << "computed: ";
      solout(sol);
      cout << "error: ";
      for (int i=0; i<s; i++) cout << fabs(sol[i]-gtsol[i]) << "  ";
      cout << endl;
      cerr << "\n\n\nthis happens if the roots of the polynomial are close together\n\n\n" 
           <<endl; 
      BIASASSERT(false);
    } 
  }
  if (debug) cout << endl << endl;
}

void errout(vector<complex<POLYNOMIALSOLVE_TYPE> >& sol, 
            vector<POLYNOMIALSOLVE_TYPE>& gtsol)
{
  sort(sol.begin(), sol.end());
  int s=(int)sol.size();
  cout.precision(20);
  if (debug) cout << "error: ";
  for (int i=0; i<s; i++){
    if (debug) cout << fabs((sol[i]).real()-gtsol[i]) << "  ";
    if (fabs((sol[i]).real()-gtsol[i])>maxerr) {
      cerr << "\n\n\nerror bigger than "<<maxerr<<"\n" <<endl; 
      debug=true;
      cout << "ground truth: ";
      solout(gtsol);
      cout << "computed: ";
      solout(sol);
      cout << "error: ";
      for (int i=0; i<s; i++) cout << fabs((sol[i]).real()-gtsol[i]) << "  ";
      cout << endl;
      cerr << "\n\n\nthis happens if the roots of the polynomial are close together\n\n\n" 
           <<endl; 
      BIASASSERT(false);
    } 
  }
  if (debug) cout << endl << endl;
}

/** @author woelk 08/2004 */
//int main(int argc, char *argv[])
int main()
{
  PolynomialSolve ps;

  //vector<double> coeff(5);
// //   coeff[0]=-3.26315441046362053612028830685;
// //   coeff[1]=13.0504510585815651779739710037;
// //   coeff[2]=-19.574778123864497558770381147;
// //   coeff[3]=13.0508215983482731559206513339;
// //   coeff[4]=-3.26334012093903513829218354658;
//   
//   coeff[0]=-5.3812376413659608331840900064;
//   coeff[1]=21.5206398747126961268349987222;
//   coeff[2]=-32.2795617934279306382450158708;
//   coeff[3]=21.5221584472618552297262795037;
//   coeff[4]=-5.38199880981732103890635698917;
//   
//   vector<double> roots;
//   ps.Quartic(coeff, roots);
//   //ps.Numeric(coeff, roots);
//   
//   /*
//   octave:2> x=roots(coeff)
//   x =
// 
//   0.99948 + 0.02686i
//   0.99948 - 0.02686i
//   1.00105 + 0.00000i
//   0.99933 + 0.00000i
//   
//   coefs2:
//   0.999231378314545 + 0.030940956057273i
//   0.999231378314545 - 0.030940956057273i
//   1.004259532432780 + 0.000000000000000i
//   0.996476651583676 + 0.000000000000000i
//   
//   */
//   
//   if(roots.size()==0)
//   {
//     cout<<"PolynomialSolve found no roots! although"<<endl;
//   }
//   else
//   {
//     for(unsigned int r=0; r<roots.size(); r++)
//     {
//     cout<<r<<":  "<<roots[r]<<endl;
//     }
//   }
  
  Random ran;

  int num=9000;
  //num=1; debug=true; ps.AddDebugLevel(D_PS_COEFF);
  int maxdeg=8;
  double maxsol=5e1;
  
  int deg, nums;
  vector<POLYNOMIALSOLVE_TYPE> coeff, sol, gtsol;
  vector<complex<POLYNOMIALSOLVE_TYPE> > csol;

  for (int n=0; n<num; n++){
    deg=ran.GetUniformDistributedInt(0, maxdeg);
    //    deg=5;
    if (debug){
      cout << "-----------------"<< setw(3) << n <<"------------------"<<endl;
      cout << "using a polynom of degree "<<deg<<endl;
    } else {
      if (n%10==0){
        cout << n << flush;
      } else {
        cout << "." << flush;
      }
    }

    gtsol.resize(deg);
    
    // first generate random solutions
    for (int d=0; d<deg; d++){
      gtsol[d]=ran.GetUniformDistributed(-maxsol, maxsol);
      //gtsol[d]=d+1;
      if (debug)
        cout << d << "th root: "<<gtsol[d]<<endl;
    }
    sort(gtsol.begin(), gtsol.end());

    // now calculate the coefficients
    ps.CalCoefficients(gtsol, coeff);
    if (debug){
      for (int i=0; i<deg+1; i++){
        cout << i << "th coefficient: "<< coeff[i]<<endl;
      }
    }
    /* // test for CheckCoefficients
    coeff[deg]=0.0;     if (deg>1) coeff[deg-1]=1e-12;
    ps.CheckCoefficients(coeff, 1e-11);
    if (debug){
      for (int i=0; i<coeff.size(); i++){
        cout << i << "th coefficient: "<< coeff[i]<<endl;
      }
    }
    */

    // check if coefficients are correct
    if (debug) {
      cout << "polynomial at ground truth roots: "<<endl;
      for (int i=0; i<deg; i++){
        cout << "p("<<-gtsol[i]<<") = "<<ps.EvaluatePolynomial(gtsol[i], coeff)
             << "   ";
        cout << endl;
      }
    }
    
    nums=ps.Solve(coeff, sol);
    if (debug) cout <<"Solve (real)    found "<<nums<<" solutions:\n";
    if (nums>0) { solout(sol); errout(sol, gtsol); }
    nums=ps.Solve(coeff, csol);
    if (debug) cout <<"Solve (complex) found "<<nums<<" solutions:\n";
    if (nums>0) { solout(csol); errout(csol, gtsol); }

    nums=ps.Numeric(coeff, sol);
    if (debug)cout<<"Numeric (real)    found "<<nums<<" solutions:\n";
    if (nums>0) { solout(sol); errout(sol, gtsol); }
    nums=ps.Numeric(coeff,csol);
    if (debug) cout<<"Numeric (complex) found "<<nums<<" solutions:\n";
    if (nums>0) { solout(csol); errout(csol, gtsol); }

    if (deg<=4){
      nums=ps.Analytic(coeff, sol);
      if (debug) cout << "Analytic (real)    found "<<nums<<" solutions:\n";
      if (nums>0) { solout(sol); errout(sol, gtsol); }
    }
    if (deg<=3){
      nums=ps.Analytic(coeff, csol);
      if (debug) cout << "Analytic (complex) found "<<nums<<" solutions:\n";
      if (nums>0) { solout(csol); errout(csol, gtsol); }
    }

  }

  return 0;
}