PolynomialSolve.hh

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


#ifndef _BIAS_POLYNOMIALSOLVE_
#define _BIAS_POLYNOMIALSOLVE_

#include <Base/Common/BIASpragmaStart.hh>

#include <cfloat>
#include <vector>
#include <complex>
#include "Base/Debug/Debug.hh"
#include "Base/Debug/Error.hh"
 
#define POLYNOMIALSOLVE_TYPE double

#define D_PS_GSL      0x00000001
#define D_PS_NLREFINE 0x00000002
#define D_PS_COEFF    0x00000004

namespace BIAS {

  /** @class PolynomialSolve
      @ingroup g_mathalgo
      @brief base class for solving polynomial equations 
      
      @author koeser 10/2003, woelk 08/2004  */
  class BIASMathAlgo_EXPORT PolynomialSolve : public Debug 
  {
    
  public:
    PolynomialSolve() {PowellCoeff_=NULL; PowellCoeffNum_=0;};
    ~PolynomialSolve() {};

    /** @brief solve polynomial of arbitrary order n
        coeff[n]x^n+...+coeff[2]x^2+coeff[1]x+coeff[0]=0
        coeff[n]!=0.0 is asserted
        Uses analytical solution if n<=4
        @return number of real solutions, negative value on error 
        @author koeser, woelk 08/2004 untested */
    int Solve(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
              std::vector<POLYNOMIALSOLVE_TYPE>& sol);

    /** @brief solve polynomial of arbitrary order n
        coeff[n]x^n+...+coeff[2]x^2+coeff[1]x+coeff[0]=0
        coeff[n]!=0.0 is asserted.
        Uses analytical solution if n<=3
        @return number of solutions, negative value on error 
        @author woelk 08/2004 untested */
    int Solve(const std::vector<POLYNOMIALSOLVE_TYPE>&  coeff,
              std::vector<std::complex<POLYNOMIALSOLVE_TYPE> >& sol);

    /** @brief solves polynomial of arbitrary order n numerically
        coeff[n]x^n+...+coeff[2]x^2+coeff[1]x+coeff[0]=0
        coeff[n]!=0.0 is asserted.
        @return number of real solutions, negative value on error 
        @author woelk 08/2004 untested */
    int Numeric(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                std::vector<POLYNOMIALSOLVE_TYPE>& sol);
    /** @brief solves polynomial of arbitrary order n numerically
        coeff[n]x^n+...+coeff[2]x^2+coeff[1]x+coeff[0]=0
        coeff[n]!=0.0 is asserted.
        @return number of complex solutions, negative value on error 
        @author woelk 08/2004 untested */
    int Numeric(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                std::vector<std::complex<POLYNOMIALSOLVE_TYPE> >& sol);

    /** @brief solves polynomial of arbitrary order n<=4 analytically
        coeff[n]x^n+...+coeff[2]x^2+coeff[1]x+coeff[0]=0
        coeff[n]!=0.0 is asserted.
        @return number of real solutions, negative value on error 
        @author woelk 08/2004 untested */
    int Analytic(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                 std::vector<POLYNOMIALSOLVE_TYPE>& sol);
    /** @brief solves polynomial of arbitrary order n<=3 analytically
        coeff[n]x^n+...+coeff[2]x^2+coeff[1]x+coeff[0]=0
        coeff[n]!=0.0 is asserted.
        @return number of complex solutions, negative value on error 
        @author woelk 08/2004 untested */
    int Analytic(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                 std::vector<std::complex<POLYNOMIALSOLVE_TYPE> >& sol);

    /** @brief solve a polynomial of degree 1
        coeff[1]*x+coeff[0]=0
        coeff[1] != 0.0 is asserted
        returns the number of real solutions (always 1)
        @author koeser, woelk 08/2004 */
    int Linear(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
               std::vector<POLYNOMIALSOLVE_TYPE>& sol);
    /** @brief solve a polynomial of degree 1
        coeff[1]*x+coeff[0]=0
        coeff[1] != 0.0 is asserted
        returns the number of complex solutions (always 1)
        @author woelk 08/2004 */
    int Linear(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
               std::vector<std::complex<POLYNOMIALSOLVE_TYPE> >& sol);

    /** @brief solve a polynomial of degree 2
        coeff[2]*x^2+coeff[1]*x+coeff[0]=0
        coeff[2] != 0.0 is asserted
        returns the number of real solutions (0, 1 or 2)
        @author koeser, woelk 08/2004         
    */
    int Quadratic(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                  std::vector<POLYNOMIALSOLVE_TYPE>& sol);
    /** @brief solve a polynomial of degree 2
        coeff[2]*x^2+coeff[1]*x+coeff[0]=0
        coeff[2] != 0.0 is asserted
        returns the number of complex solutions (alway 2)
        @author woelk 08/2004 */
    int Quadratic(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                  std::vector<std::complex<POLYNOMIALSOLVE_TYPE> >& sol);

    /** @brief solve a polynomial of degree 3
        coeff[3]*x^3+coeff[2]*x^2+coeff[1]*x+coeff[0]=0
         coeff[3] != 0.0 is asserted
         returns the number of real solutions 
         @author koeser, woelk 08/2004 */
    int Cubic(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
              std::vector<POLYNOMIALSOLVE_TYPE>& sol);
    /** @brief solve a polynomial of degree 3
        coeff[3]*x^3+coeff[2]*x^2+coeff[1]*x+coeff[0]=0
        coeff[3] != 0.0 is asserted
        returns the number of real solutions 
        @author woelk 08/2004 */
    int Cubic(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
              std::vector<std::complex<POLYNOMIALSOLVE_TYPE> >& sol);

    /** @brief solve a polynomial of degree 4
        coeff[4]*x^4+coeff[3]*x^3+coeff[2]*x^2+coeff[1]*x+coeff[0]=0
        coeff[4] != 0.0 is asserted
        returns the number of real solutions 
        @author koeser, woelk 08/2004 
    ATTENTION: under certain circumstances Quartic() doesn't seem to work
               correctly and does not return any solutions.
    EXAMPLE  : for the following polynom   
          coeff[0]=-5.3812376413659608331840900064;
          coeff[1]=21.5206398747126961268349987222;
          coeff[2]=-32.2795617934279306382450158708;
          coeff[3]=21.5221584472618552297262795037;
          coeff[4]=-5.38199880981732103890635698917;
    
   Quartic() does not return any roots. 
        Expected roots are:
            0.999231378314545 + 0.030940956057273i
            0.999231378314545 - 0.030940956057273i
            1.004259532432780 + 0.000000000000000i
            0.996476651583676 + 0.000000000000000i
    */
    int Quartic(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                std::vector<POLYNOMIALSOLVE_TYPE>& sol);

    /** @brief removes leading zeros in coefficients
        @author woelk 08/2004 */
    void CheckCoefficients(std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                           double eps=DBL_EPSILON);


    /** @brief inverse function, calculates coefficients from the solutions 
        @author woelk 08/2004 */
    void CalCoefficients(const std::vector<POLYNOMIALSOLVE_TYPE>& sol,
                         std::vector<POLYNOMIALSOLVE_TYPE>& coeff);

    /** @brief numerically robust way to evaluate a polynomial at position x,
        uses Horner scheme (same principle as in gsl_poly_eval)
        @author koeser */
    inline POLYNOMIALSOLVE_TYPE 
    EvaluatePolynomial(const POLYNOMIALSOLVE_TYPE x, 
                       const std::vector<POLYNOMIALSOLVE_TYPE>& coeff) const;

    /** @brief determine number of real solutions of polynomial */
    int GetNumberOfRealSolutions(const std::vector<POLYNOMIALSOLVE_TYPE>& 
                                 coeff);

    /** @brief determine whether the polynomial has a solution in IR */
    bool HasRealSolution(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff)
    { return (GetNumberOfRealSolutions(coeff)>0); }


    /** @brief does a non linear refinement using Powel from minpack
        for real roots (non imaginary) only 
        @author woelk, koeser 08/2004 */
    int NonLinearRefine(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                        std::vector<std::complex<POLYNOMIALSOLVE_TYPE> >& sol);

    /** @brief does a non linear refinement using Powel from minpack
        for real roots (non imaginary) only 
        @author woelk, koeser 11/2007 */
    int NonLinearRefine(const std::vector<POLYNOMIALSOLVE_TYPE>& coeff,
                        std::vector<POLYNOMIALSOLVE_TYPE>& sol);

    /** @brief given locations x and measurements y=f(x) approximate f by a 
     *  polynomial of arbitrary degree and return coefficients
     *
     *  least squares approximation of parameters a_i of a polynomial
     *  which minimizes (a_0 + a_1x + a_2x^2 + ... + a_degree x^degree - y)^2
     *  summed over all pairs (x;y), where i runs from 0 to degree
     *  @return 0 on success, <0 for degenerate situation, >0 optimizer error
     *  @author koeser 04/2007 */
    int FitPolynomial(const unsigned int degree,
                      const std::vector<double>& x,
                      const std::vector<double>& y,
                      std::vector<double> &coefficients);

  public:       
    int PowellCoeffNum_;
    double *PowellCoeff_;

  protected:
    /** helper function for CalCoefficients
        @author woelk 08/2004 */
    POLYNOMIALSOLVE_TYPE 
    _GetCoeff(int order, int b, const std::vector<POLYNOMIALSOLVE_TYPE>& sol);
 
  };

  //////////////////////////////////////////////////////////////////////
  //        inline implementation
  //////////////////////////////////////////////////////////////////////

  inline POLYNOMIALSOLVE_TYPE 
  PolynomialSolve::EvaluatePolynomial(const POLYNOMIALSOLVE_TYPE x, const
                                      std::vector<POLYNOMIALSOLVE_TYPE>& coeff) const
  {
    register int n = (int)coeff.size() - 1;
    POLYNOMIALSOLVE_TYPE result = coeff[n--];
    while (n>=0) result = result*x + coeff[n--];
    return result; 
  }


} // end namespace

#include <Base/Common/BIASpragmaEnd.hh>

#endif