ossimPolynom< T, DIM > Class Template Reference

template class for multivariate polynom algebra More...

#include <ossimPolynom.h>

Classes
struct	EXP_TUPLE_LESSTHAN
	type to store multivariate input More...

Public Types
typedef std::vector< int >	EXP_TUPLE
	inner types More...
typedef std::vector< T >	VAR_TUPLE
	type to store exponent tuples More...
typedef std::map< EXP_TUPLE, T, EXP_TUPLE_LESSTHAN >	MONOM_MAP
typedef std::set< EXP_TUPLE, EXP_TUPLE_LESSTHAN >	EXPT_SET
	for storing polynom More...

Public Member Functions
	ossimPolynom (const T &epsilon=0)
	for storing set of exponent tuples More...
	ossimPolynom (const ossimPolynom &p)
	~ossimPolynom ()
const ossimPolynom &	operator= (const ossimPolynom< T, DIM > &pt)
void	setMonom (const EXP_TUPLE &m, const T &v)
void	setMonom (const int mexp[DIM], const T &v)
void	delMonom (const EXP_TUPLE &m)
T	getCoeff (const EXP_TUPLE &m) const
void	nullify ()
bool	isNull () const
const MONOM_MAP &	getMonoms () const
const T &	getEpsilon () const
	no setEpsilon beacause might erase monoms More...
bool	operator== (const ossimPolynom &pt) const
	comparison operators -don't compare theEpsilon values -use my own epsilon in comparisons (not the compared to's epsilon) More...
bool	operator!= (const ossimPolynom &pt) const
bool	isNegligible (const T &v) const
	can v be considered as zero? More...
int	getOrder (int d) const
	orders More...
int	getTotalOrder () const
	returns maximum order of monoms More...
T	eval (const VAR_TUPLE &v) const
	evaluation : needs DIM values as input More...
void	pdiff (int pdim, ossimPolynom &result) const
	partial differentiation polynom More...
const ossimPolynom &	operator*= (const T &k)
	operations with scalar More...
ossimPolynom	operator+ (const ossimPolynom &p) const
	arithmetic operators More...
ossimPolynom	operator- (const ossimPolynom &p) const
const ossimPolynom &	operator+= (const ossimPolynom &p)
const ossimPolynom &	operator-= (const ossimPolynom &p)
ossimPolynom	operator* (const ossimPolynom &p) const
	product operator : use epsilon from left side More...
const ossimPolynom &	operator*= (const ossimPolynom &p)
std::ostream &	print (std::ostream &os) const
	I/O. More...
std::ostream &	printNice (std::ostream &os, const char symbols[DIM])
	classic representation (bad accuracy, for display only) More...
std::istream &	import (std::istream &is)
	note that it can only import for the template type T and dimesnion DIM More...
EXPT_SET	builExpSet (const EXP_TUPLE &orders) const
	constructs simple exponent tuples set for using LMSfit need order for each dimension More...
void	addExpTupleRight (int newdim, int totalOrder, EXPT_SET &eset) const
	concatenate exponents (at the right) to existing tuple set, for a given maximum total order eg: with eset={(0,1),(0,0)} , then addExpTuple(2,1,eset) = {(0,1,0,0),(0,1,0,1),(0,1,1,0), (0,0,0,0),(0,0,0,1),(0,0,1,0)} More...
bool	LMSfit (const EXPT_SET &expset, const std::vector< VAR_TUPLE > obs_input, const std::vector< T > obs_output, T *prms=NULL)
	fits the polynom to the observations using Least Mean Squares returns true on success (can fail if not enough observations) More...
bool	SLSfit (const EXPT_SET &expset, const std::vector< VAR_TUPLE > obs_input, const std::vector< T > obs_output, T *prms=NULL)
	Standard least squares Modified version of LMSfit that uses standard NEWMAT inverse as alternative to SVD solution. More...

Static Public Member Functions
static void	addExpTuple (EXP_TUPLE &m1, const EXP_TUPLE &m2)

Protected Member Functions
void	discardNullMonoms ()
	positive values below epsilon are considered 0 More...
void	addMonomQuick (const EXP_TUPLE &m, const T &v)
	add value without testing for negligible More...
NEWMAT::Matrix	invert (const NEWMAT::Matrix &m) const
	invert stolen from ossimRpcSolver More...

Protected Attributes
MONOM_MAP	theMonoms
	protected data members More...
T	theEpsilon
	associate a scalar to each tuple of orders : monom More...

Detailed Description

template<class T, int DIM = 1>
class ossimPolynom< T, DIM >

template class for multivariate polynom algebra

T : the storage type, constraints: must have 0 (zero) as value, support ops fabs + - * / DIM : the dimension of the input space, integer (>=1, default : 1)

stores a set of monoms, a monom is (exponent tuples + coefficient) requires a precion (epsilon) for comparisons note: monoms absolute values below epsilon are removed from the map

Definition at line 38 of file ossimPolynom.h.

Member Typedef Documentation

EXP_TUPLE

template<class T, int DIM = 1>

typedef std::vector< int > ossimPolynom< T, DIM >::EXP_TUPLE

inner types

Definition at line 44 of file ossimPolynom.h.

EXPT_SET

template<class T, int DIM = 1>

typedef std::set< EXP_TUPLE, EXP_TUPLE_LESSTHAN > ossimPolynom< T, DIM >::EXPT_SET

for storing polynom

Definition at line 68 of file ossimPolynom.h.

MONOM_MAP

template<class T, int DIM = 1>

typedef std::map< EXP_TUPLE, T , EXP_TUPLE_LESSTHAN > ossimPolynom< T, DIM >::MONOM_MAP

Definition at line 66 of file ossimPolynom.h.

VAR_TUPLE

template<class T, int DIM = 1>

typedef std::vector< T > ossimPolynom< T, DIM >::VAR_TUPLE

type to store exponent tuples

Definition at line 45 of file ossimPolynom.h.

Constructor & Destructor Documentation

ossimPolynom() [1/2]

template<class T, int DIM = 1>

ossimPolynom< T, DIM >::ossimPolynom ( const T &  epsilon = 0 )

inline

for storing set of exponent tuples

construction : must supply epsilon value, default 0

Definition at line 75 of file ossimPolynom.h.

      : theEpsilon(epsilon)
   {}

ossimPolynom() [2/2]

template<class T, int DIM = 1>

ossimPolynom< T, DIM >::ossimPolynom ( const ossimPolynom< T, DIM > &  p )

inline

Definition at line 79 of file ossimPolynom.h.

                                       : 
      theMonoms(p.getMonoms()),
      theEpsilon(p.getEpsilon())
   {}

~ossimPolynom()

template<class T, int DIM = 1>

ossimPolynom< T, DIM >::~ossimPolynom ( )

inline

Definition at line 84 of file ossimPolynom.h.

   {}

Member Function Documentation

addExpTuple()

template<class T, int DIM = 1>

static void ossimPolynom< T, DIM >::addExpTuple ( EXP_TUPLE &  m1, const EXP_TUPLE &  m2  )

inlinestatic

Definition at line 344 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::operator*().

   {
      for(int i=0;i<DIM;++i) {
         m1[i] += m2[i];
      }
   }

addExpTupleRight()

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::addExpTupleRight ( int  newdim, int  totalOrder, EXPT_SET &  eset  ) const

inline

concatenate exponents (at the right) to existing tuple set, for a given maximum total order eg: with eset={(0,1),(0,0)} , then addExpTuple(2,1,eset) = {(0,1,0,0),(0,1,0,1),(0,1,1,0), (0,0,0,0),(0,0,0,1),(0,0,1,0)}

Definition at line 569 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::addExpTupleRight(), and ossimPolynomProjection::setupDesiredExponents().

 {
    EXPT_SET newset;
    // add a copy off eset for each order with the specific last dim
    for(int o=0; o<=totalOrder ; ++o)
    {
       EXPT_SET extset; //extended set
       if (eset.size()==0)
       {
          EXP_TUPLE tu(1);
          tu[0]=o;
          extset.insert(tu);
       } else {
          //we have to construct a new set from eset, extending dimension
          // cause: stored tuples cannot be compared at different dimensions
          for(typename EXPT_SET::iterator sit = eset.begin(); sit != eset.end(); ++sit)
          {
             EXP_TUPLE tu(*sit);
             tu.push_back(o);
             extset.insert(tu);
          }
       }
       //recursively add remaining dimensions
       if (newdim>1)
       {
          addExpTupleRight(newdim-1, totalOrder-o, extset); //only dimension decreases
       }
       //add full set for the specific order
       newset.insert(extset.begin(), extset.end());
    }
    eset=newset; //overwrite
 }

addMonomQuick()

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::addMonomQuick ( const EXP_TUPLE &  m, const T &  v  )

inlineprotected

add value without testing for negligible

Definition at line 807 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::operator*().

   {
      typename MONOM_MAP::iterator it = theMonoms.find(m);
      if (it != theMonoms.end())
      {
         it->second  += v;
      } else {
         theMonoms.insert( MONOM_MAP::value_type(m,v));
      }
   }

builExpSet()

template<class T, int DIM = 1>

EXPT_SET ossimPolynom< T, DIM >::builExpSet ( const EXP_TUPLE &  orders ) const

inline

constructs simple exponent tuples set for using LMSfit need order for each dimension

Definition at line 534 of file ossimPolynom.h.

 {
    EXPT_SET eset;
    if (orders.size() != DIM)
    {
       std::cerr<<"ossimPolynom::import ERROR bad dimension for parameter, need "<<DIM<<" elements"<<std::endl;
       return eset;
    }
    //initialise variable exponent tuple
    EXP_TUPLE et(DIM);
    for(int d=0;d<DIM;++d) et[d]=0;
    while (et[0] <= orders[0])
    {
      //add tuple to set
      eset.insert(et);

      //increment tuple within bounds
      et[DIM-1]++;
      for(int d=DIM-1 ; d>=0 ; --d)
      {
         if ((et[d] > orders[d]) && (d>0))
         {
            et[d] = 0;
            et[d-1]++;
         }
      }
    }
    return eset;
 }

delMonom()

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::delMonom ( const EXP_TUPLE &  m )

inline

Definition at line 119 of file ossimPolynom.h.

   {
      theMonoms.erase(m); //TBC TBD: what happens if m is not in the map?
   }

discardNullMonoms()

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::discardNullMonoms ( )

inlineprotected

positive values below epsilon are considered 0

method to erase all negligible monoms : user don't need that (automatic)

Definition at line 789 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::operator*(), and ossimPolynom< ossim_float64, 3 >::operator*=().

   {     
      std::vector< typename MONOM_MAP::iterator > erasev; //storage for iterators on elements to erase
      for (typename MONOM_MAP::iterator it = theMonoms.begin(); it != theMonoms.end() ; ++it )
      {
         if (isNegligible(it->second)) erasev.push_back(it);
      }
      //erase all elements afterwards
      for ( typename std::vector< typename MONOM_MAP::iterator >::const_iterator vit = erasev.begin(); vit != erasev.end() ; ++vit )
      {
         theMonoms.erase(*vit); //*vit is  an iterator in theMonoms
      }      
   }

eval()

template<class T, int DIM = 1>

T ossimPolynom< T, DIM >::eval ( const VAR_TUPLE &  v ) const

inline

evaluation : needs DIM values as input

Definition at line 219 of file ossimPolynom.h.

Referenced by ossimPolynomProjection::lineSampleHeightToWorld(), and ossimPolynomProjection::worldToLineSample().

   {
      T ev = 0;
      //loop on monoms. TBD optimize powers using map order
      typename MONOM_MAP::const_iterator it;
      int p;
      for ( it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {
         //compute powers
         T mv = it->second;
         for(int d=0;d<DIM;++d)
         {
            p = (it->first)[d];
            if (p != 0)
            {
               mv *= std::pow( v[d], p );
            }
         }
         //add momom value
         ev += mv;
      }
      return ev;
   }

getCoeff()

template<class T, int DIM = 1>

T ossimPolynom< T, DIM >::getCoeff ( const EXP_TUPLE &  m ) const

inline

Definition at line 124 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::operator+=(), ossimPolynom< ossim_float64, 3 >::operator-=(), and ossimPolynom< ossim_float64, 3 >::operator==().

   {
      typename MONOM_MAP::const_iterator it = theMonoms.find(m);
      if (it != theMonoms.end())
      {
         return it->second;
      } else {
         return 0;
      }
   }

getEpsilon()

template<class T, int DIM = 1>

const T& ossimPolynom< T, DIM >::getEpsilon ( ) const

inline

no setEpsilon beacause might erase monoms

Definition at line 150 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::invert(), ossimPolynom< ossim_float64, 3 >::operator*(), ossimPolynom< ossim_float64, 3 >::operator=(), and ossimPolynom< ossim_float64, 3 >::print().

   {
      return theEpsilon;
   }

getMonoms()

template<class T, int DIM = 1>

const MONOM_MAP& ossimPolynom< T, DIM >::getMonoms ( ) const

inline

Definition at line 145 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::eval(), ossimPolynom< ossim_float64, 3 >::getOrder(), ossimPolynom< ossim_float64, 3 >::getTotalOrder(), ossimPolynom< ossim_float64, 3 >::import(), ossimPolynom< ossim_float64, 3 >::operator*(), ossimPolynom< ossim_float64, 3 >::operator*=(), ossimPolynom< ossim_float64, 3 >::operator+=(), ossimPolynom< ossim_float64, 3 >::operator-=(), ossimPolynom< ossim_float64, 3 >::operator=(), ossimPolynom< ossim_float64, 3 >::operator==(), ossimPolynom< ossim_float64, 3 >::pdiff(), ossimPolynom< ossim_float64, 3 >::print(), and ossimPolynom< ossim_float64, 3 >::printNice().

   {
      return theMonoms;
   }

getOrder()

template<class T, int DIM = 1>

int ossimPolynom< T, DIM >::getOrder ( int  d ) const

inline

orders

returns maximum order of monoms for a specific dimension (d starts at 0)

Definition at line 186 of file ossimPolynom.h.

   {
      if ((d>=DIM) || (d<0)) return -1; //error = no dimension

      // loop on monoms
      int order = -1; //for null polynom
      int corder;
      typename MONOM_MAP::const_iterator it;
      for ( it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {
         corder = (it->first)[d];
         if ( corder > order ) order = corder;
      }
      return order;
   }

getTotalOrder()

template<class T, int DIM = 1>

int ossimPolynom< T, DIM >::getTotalOrder ( ) const

inline

returns maximum order of monoms

Definition at line 202 of file ossimPolynom.h.

   {
      int order = -1; //for null polynom
      int sorder;
      typename MONOM_MAP::const_iterator it;
      for ( it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {
         sorder = 0;
         for (int d=0;d<DIM;++d) sorder+=(it->first)[d];
         if ( sorder > order ) order = sorder;
      }
      return order;
   }

import()

template<class T, int DIM = 1>

std::istream& ossimPolynom< T, DIM >::import ( std::istream &  is )

inline

note that it can only import for the template type T and dimesnion DIM

Definition at line 440 of file ossimPolynom.h.

Referenced by operator>>().

   {
      nullify();

      //extract data in brackets [ ]
      //swallow bracket
      ossimString tempString;
      char tempChar;
      is.get(tempChar);
      if (!is) 
      {
         std::cerr<<"ossimPolynom::import ERROR, cannot read left bracket ["<<std::endl;
         return is;
      }
      if (tempChar != '[')
      {
         std::cerr<<"ossimPolynom::import ERROR, missing left bracket ["<<std::endl;
         return is;
      }
      const int BUFSIZE=32768;  //should be enough fro most apps (TBC TBD : allow loops if not enough)
      char buffer[BUFSIZE];
      is.getline(buffer, BUFSIZE, ']');
      if (!is)
      {
         std::cerr<<"ossimPolynom::import ERROR, cannot read after left bracket ["<<std::endl;
         return is;
      }     
      if (is.gcount() >= BUFSIZE-1)
      {
         std::cerr<<"ossimPolynom::import ERROR, cannot find right bracket ] after "<<BUFSIZE-1<<"characters"<<std::endl;
         return is;
      }
      tempString = buffer; //no more brackets

      //split string using semicolons
      std::vector< ossimString > subparts = tempString.explode(";");
      //loop on subparts
      for (typename std::vector< ossimString >::const_iterator it=subparts.begin() ; it!=subparts.end() ; ++it )
      {
         ossimString sp = it->trim();
         //check for epsilon
         ossimString aft_eps=sp.after("eps=");
         if (aft_eps.size()>0)
         {
            //get epsilon value
            aft_eps.trim();
            theEpsilon = static_cast< T >(aft_eps.toDouble());
         } else {
            //no epsilon : must be a monom
            ossimString scalpart=sp.before("(");
            if (scalpart.size() < sp.size())
            {
              T coeff = static_cast< T >(scalpart.toDouble());
              ossimString expopart = ((sp.before(")")).after("(")).trim();
              if (expopart.size() == 0)
              {
                 std::cerr<<"ossimPolynom::import ERROR, cannot find ) in monom or empty monom"<<std::endl;
                 return is;
              }
              std::vector< ossimString > vexpo = expopart.explode(",");
              if (vexpo.size() != DIM)
              {
                 std::cerr<<"ossimPolynom::import ERROR, bad number of exponents in monom (need "<<DIM<<"): "<<vexpo.size()<<std::endl;
                 return is;
              }
              //store all exponents
              EXP_TUPLE expt(DIM);
              int d;
              std::vector< ossimString >::const_iterator eit;
              for (eit=vexpo.begin() , d=0 ; eit != vexpo.end() ; ++eit, ++d )
              {
                 expt[d] = eit->toInt(); //TBD : could check that value is integer, but how?
              }
              //add new monom (if duplicate...error)
              if (getMonoms().find(expt) != getMonoms().end())
              {
                 std::cerr<<"ossimPolynom::import ERROR, duplicate exponent tuple in polynom"<<std::endl;
                 return is;
              }
              theMonoms[expt] = coeff;

            } else {
               std::cerr<<"ossimPolynom::import ERROR, cannot find left parenthesis ( in monom "<<std::endl;
               return is;
            }
         }
      }

      return is;
   }

invert()

template<class T, int DIM = 1>

NEWMAT::Matrix ossimPolynom< T, DIM >::invert ( const NEWMAT::Matrix &  m ) const

inlineprotected

invert stolen from ossimRpcSolver

Definition at line 822 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::LMSfit().

   {
      ossim_uint32 idx = 0;
      NEWMAT::DiagonalMatrix d;
      NEWMAT::Matrix u;
      NEWMAT::Matrix v;

      // decompose m.t*m which is stored in Temp into the singular values and vectors.
      //
      NEWMAT::SVD(m, d, u, v, true, true);
      
      // invert the diagonal
      // this is just doing the reciprical fo all diagonal components and store back int
      // d.  ths compute d inverse.
      //
      for(idx=0; idx < (ossim_uint32)d.Ncols(); ++idx)
      {
         if(d[idx] > getEpsilon()) //adpated here for epsilon
         {
            d[idx] = 1.0/d[idx];
         }
         else
         {
            d[idx] = 0.0;
         }
      }

      //compute inverse of decomposed m;
      return v*d*u.t();
   }

isNegligible()

template<class T, int DIM = 1>

bool ossimPolynom< T, DIM >::isNegligible ( const T &  v ) const

inline

can v be considered as zero?

Definition at line 178 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::discardNullMonoms(), ossimPolynom< ossim_float64, 3 >::operator==(), and ossimPolynom< ossim_float64, 3 >::setMonom().

   {
      return ( fabs(v) <= theEpsilon );
   }   

isNull()

template<class T, int DIM = 1>

bool ossimPolynom< T, DIM >::isNull ( ) const

inline

Definition at line 140 of file ossimPolynom.h.

   {
      return (theMonoms.size() == 0);
   }

LMSfit()

template<class T, int DIM = 1>

bool ossimPolynom< T, DIM >::LMSfit ( const EXPT_SET &  expset, const std::vector< VAR_TUPLE >  obs_input, const std::vector< T >  obs_output, T *  prms = NULL  )

inline

fits the polynom to the observations using Least Mean Squares returns true on success (can fail if not enough observations)

• also updates rms error(root mean square) NOTES: inputs must have same size and must be ordered the same way use builExpSet() to construct classic polynoms TODO: add weights to observations

Definition at line 611 of file ossimPolynom.h.

Referenced by ossimPolynomProjection::optimizeFit().

{
   //init
   nullify();

   //check size
   int nobs = (int)obs_input.size();
   if (nobs != (int)obs_output.size())
   {
      std::cerr<<"ossimPolynom::LMSfit ERROR observation input/output must have the same size"<<std::endl;
      return false;
   }
   if (nobs<=0)
   {
      std::cerr<<"ossimPolynom::LMSfit ERROR observation count is zero"<<std::endl;
      return false;
   }
   int ncoeff = (int)expset.size();  
   if (ncoeff<=0)
   {
      std::cerr<<"ossimPolynom::LMSfit ERROR exponent count is zero"<<std::endl;
      return false;
   }
   
   //construct LMS linear system (using OSSIM matrices)
   // M.k = v
   // M : monom matrix
   // k : polynbm coefficients
   // v : output_obs
   NEWMAT::Matrix M(nobs, ncoeff);
   double elt;
   int p;
   typename EXPT_SET::const_iterator cit;
   typename std::vector< VAR_TUPLE >::const_iterator oit;
   int i,j;
   for (cit=expset.begin(), j=1; cit != expset.end() ; ++cit, ++j) //loop on exponent tuples
   {
      for(oit=obs_input.begin(), i=1; oit!=obs_input.end();++oit, ++i) //loop on observations
      {
         //compute powers using observation position
         elt=1.0;
         for(int d=0;d<DIM;++d)
         {
            p = (*cit)[d];
            if (p != 0)
            {
               elt *= std::pow( (*oit)[d], p );
            }
         }
         //init M
         M(i,j) = elt; //NEWMAT indices start at 1
      }
   }

   NEWMAT::ColumnVector v(nobs);
   for(int o=0;o<nobs;++o)
   {
      v(o+1) = obs_output[o];
   }

   //build LMS symmetric matrix tM.M
   //build best fit
   NEWMAT::ColumnVector bfit = invert(M.t()*M)*M.t()*v; //TBD : check inversion

   //compute RMS (optional, if rms non null)
   if (prms!=NULL)
   {
      NEWMAT::ColumnVector delta = M*bfit - v;
      *prms = std::sqrt( delta.SumSquare() / nobs);
   }

   //init polynom
   for (cit=expset.begin(), j=1; cit != expset.end() ; ++cit, ++j) //loop on exponent tuples
   {
      setMonom(*cit, bfit(j));
   }
   return true;
}

nullify()

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::nullify ( )

inline

Definition at line 135 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::import(), ossimPolynom< ossim_float64, 3 >::LMSfit(), ossimPolynom< ossim_float64, 3 >::pdiff(), and ossimPolynom< ossim_float64, 3 >::SLSfit().

   {
      theMonoms.clear();
   }

operator!=()

template<class T, int DIM = 1>

bool ossimPolynom< T, DIM >::operator!= ( const ossimPolynom< T, DIM > &  pt ) const

inline

Definition at line 173 of file ossimPolynom.h.

   {
      return !operator==(pt);
   }

operator*()

template<class T, int DIM = 1>

ossimPolynom ossimPolynom< T, DIM >::operator* ( const ossimPolynom< T, DIM > &  p ) const

inline

product operator : use epsilon from left side

Definition at line 315 of file ossimPolynom.h.

   {
      //do a stupid multiplication (sum all monom pairs)      
      ossimPolynom< T , DIM > prod(getEpsilon());
      T coeff;
      //loop on p monoms
      for ( typename MONOM_MAP::const_iterator it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {      
         for ( typename MONOM_MAP::const_iterator pit = p.getMonoms().begin(); pit != p.getMonoms().end() ; ++pit )
         {
            coeff = it->second * pit->second;
            if (coeff!=0)
            {
               EXP_TUPLE prodExp(it->first);
               addExpTuple(prodExp, pit->first);
               prod.addMonomQuick(prodExp, coeff);
            }
         }
      }
      //scan for null monoms and discard
      prod.discardNullMonoms();
      return prod;
   }

operator*=() [1/2]

template<class T, int DIM = 1>

const ossimPolynom& ossimPolynom< T, DIM >::operator*= ( const T &  k )

inline

operations with scalar

Definition at line 266 of file ossimPolynom.h.

   {
      //loop on monoms
      for ( typename MONOM_MAP::const_iterator it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {
         it->second *= k;
      }
      discardNullMonoms();
   }

operator*=() [2/2]

template<class T, int DIM = 1>

const ossimPolynom& ossimPolynom< T, DIM >::operator*= ( const ossimPolynom< T, DIM > &  p )

inline

Definition at line 339 of file ossimPolynom.h.

   {
      return operator=( this->operator*(p) );
   }

operator+()

template<class T, int DIM = 1>

ossimPolynom ossimPolynom< T, DIM >::operator+ ( const ossimPolynom< T, DIM > &  p ) const

inline

arithmetic operators

Definition at line 279 of file ossimPolynom.h.

   {
      ossimPolynom< T , DIM > sum(*this);
      return (sum+=p);
   }

operator+=()

template<class T, int DIM = 1>

const ossimPolynom& ossimPolynom< T, DIM >::operator+= ( const ossimPolynom< T, DIM > &  p )

inline

Definition at line 290 of file ossimPolynom.h.

   {
      typename MONOM_MAP::const_iterator it;
      //loop on p monoms
      for ( it = p.getMonoms().begin(); it != p.getMonoms().end() ; ++it )
      {
         setMonom( it->first, getCoeff(it->first) + it->second );
      }
      return *this;
   }

operator-()

template<class T, int DIM = 1>

ossimPolynom ossimPolynom< T, DIM >::operator- ( const ossimPolynom< T, DIM > &  p ) const

inline

Definition at line 284 of file ossimPolynom.h.

   {
      ossimPolynom< T , DIM > diff(*this);
      return (diff-=p);
   }

operator-=()

template<class T, int DIM = 1>

const ossimPolynom& ossimPolynom< T, DIM >::operator-= ( const ossimPolynom< T, DIM > &  p )

inline

Definition at line 301 of file ossimPolynom.h.

   {
      typename MONOM_MAP::const_iterator it;
      //loop on p monoms
      for ( it = p.getMonoms().begin(); it != p.getMonoms().end() ; ++it )
      {
         setMonom( it->first, getCoeff(it->first) - it->second );
      }
      return *this;
   }

operator=()

template<class T, int DIM = 1>

const ossimPolynom& ossimPolynom< T, DIM >::operator= ( const ossimPolynom< T, DIM > &  pt )

inline

Definition at line 87 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::operator*=().

   {
      if (this != &pt)
      {
         theEpsilon = pt.getEpsilon();
         theMonoms  = pt.getMonoms();
      }
      return *this;
   }

operator==()

template<class T, int DIM = 1>

bool ossimPolynom< T, DIM >::operator== ( const ossimPolynom< T, DIM > &  pt ) const

inline

comparison operators -don't compare theEpsilon values -use my own epsilon in comparisons (not the compared to's epsilon)

Definition at line 160 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::operator!=().

   {
      if (getMonoms().size() != pt.getMonoms().size()) return false;

      // loop on my monoms
      typename MONOM_MAP::const_iterator it;
      for ( it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {
         if ( !isNegligible(it->second - pt.getCoeff(it->first)) ) return false;
      }
      return true; //same number of identical monoms
   }

pdiff()

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::pdiff ( int  pdim, ossimPolynom< T, DIM > &  result  ) const

inline

partial differentiation polynom

Definition at line 246 of file ossimPolynom.h.

Referenced by ossimPolynomProjection::buildLineDerivatives(), and ossimPolynomProjection::buildSampDerivatives().

   {
      result.nullify();
      int ord;
      //loop on monoms
      for ( typename MONOM_MAP::const_iterator it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {
         ord = it->first[pdim];
         if (ord>=1)
         {
            EXP_TUPLE expDiff(it->first);
            expDiff[pdim] -= 1;
            result.setMonom(expDiff, it->second * ord);
         }
      }
   }

print()

template<class T, int DIM = 1>

std::ostream& ossimPolynom< T, DIM >::print ( std::ostream &  os ) const

inline

I/O.

stream serialization format : [ k1 (e1_1,e1_2,...,e1_DIM) ; k2 (e2_1,e2_2,..,e2_DIM) ; kN (eN_1,...eN_DIM)]

N is the number of monoms ei_j is exponent for dimension j and monom i

order is not important all monoms should have the same dimension : DIM you should add eps=xxxx at the beginning, separated by semi-colon ; (by default epsilon is 0)

examples: [ ] is the null polynom, [ eps=1.0e-5 ] too [ 1.0 (0) 3.5 (1) ] is polynom 1.0 + 3.5*x, with epsilon = 0 [ eps=1.0E-12 ; 2.0 (1,1,0) ; -1.0 (0,0,1) ] is polynom 2*x*y-z, with epsilon=10^-12

Definition at line 370 of file ossimPolynom.h.

   {
      static const char* monom_sep = " ; ";
      static const char* no_sep    = "";

      const char* use_sep = no_sep;

      os<<"[";
      os<<std::setprecision(16); //16 for double, TBD TBC: adapt to epsilon
      os<<std::setiosflags(std::ios_base::scientific);
      //output epsilon if not null
      if (getEpsilon() > 0)
      {         
         os<<use_sep<<"eps="<<getEpsilon();
         use_sep=monom_sep;
      }
      //loop on monoms
      for ( typename MONOM_MAP::const_iterator it = getMonoms().begin(); it != getMonoms().end() ; ++it )
      {
         os<<use_sep<<it->second<<"(";
         for(int d=0 ; d<DIM ; ++d)
         {
            /* if (d>0) */
            /* { */
            /*    os<<","; */
            /* } */
            os<<(it->first)[d];
         }
         os<<")";
         use_sep=monom_sep;
      }
      os<<"]";
      return os;
   }

printNice()

template<class T, int DIM = 1>

std::ostream& ossimPolynom< T, DIM >::printNice ( std::ostream &  os, const char  symbols[DIM]  )

inline

classic representation (bad accuracy, for display only)

Definition at line 406 of file ossimPolynom.h.

   {
      if (getMonoms().size() == 0)
      {
         os<<"0"; //zero
      } else {
         os<<std::setiosflags(std::ios_base::fixed);
         os<<std::setprecision(14); //14 for double, TBD TBC: adapt to epsilon
         //loop on monoms (map order)
         for ( typename MONOM_MAP::const_iterator it = getMonoms().begin(); it != getMonoms().end() ; ++it )
         {
            T coeff = it->second;
            if (coeff>0)
            {
               os<<"+";
            }
            os<<coeff;
            for(int d=0;d<DIM;++d) 
            {
               int ord=(it->first)[d];
               if (ord>0)
               {
                  os<<""<<symbols[d];
                  if (ord != 1)
                  {
                     os<<ord;
                  }
               }
            }
         }
      }
      return os;
   }

setMonom() [1/2]

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::setMonom ( const EXP_TUPLE &  m, const T &  v  )

inline

Definition at line 97 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::LMSfit(), ossimPolynom< ossim_float64, 3 >::operator+=(), ossimPolynom< ossim_float64, 3 >::operator-=(), ossimPolynom< ossim_float64, 3 >::pdiff(), and ossimPolynom< ossim_float64, 3 >::SLSfit().

   {
      if (isNegligible(v))
      {
         theMonoms.erase(m); //TBC TBD: what happens if m is not in the map?
      } else {
         theMonoms[m] = v;
      }
   }

setMonom() [2/2]

template<class T, int DIM = 1>

void ossimPolynom< T, DIM >::setMonom ( const int  mexp[DIM], const T &  v  )

inline

Definition at line 107 of file ossimPolynom.h.

   {
      EXP_TUPLE mexpV(mexp,mexp+DIM);
      if (isNegligible(v))
      {
         theMonoms.erase(mexpV); //TBC TBD: what happens if m is not in the map?
      } else {
         theMonoms[mexpV] = v;
      }
   }

SLSfit()

template<class T, int DIM = 1>

bool ossimPolynom< T, DIM >::SLSfit ( const EXPT_SET &  expset, const std::vector< VAR_TUPLE >  obs_input, const std::vector< T >  obs_output, T *  prms = NULL  )

inline

Standard least squares Modified version of LMSfit that uses standard NEWMAT inverse as alternative to SVD solution.

Definition at line 700 of file ossimPolynom.h.

{
   //init
   nullify();

   //check size
   int nobs = (int)obs_input.size();
   if (nobs != (int)obs_output.size())
   {
      std::cerr<<"ossimPolynom::LMSfit ERROR observation input/output must have the same size"<<std::endl;
      return false;
   }
   if (nobs<=0)
   {
      std::cerr<<"ossimPolynom::LMSfit ERROR observation count is zero"<<std::endl;
      return false;
   }
   int ncoeff = (int)expset.size();  
   if (ncoeff<=0)
   {
      std::cerr<<"ossimPolynom::LMSfit ERROR exponent count is zero"<<std::endl;
      return false;
   }
   
   // M : monom matrix
   // k : polynomial coefficients
   // v : output_obs
   NEWMAT::Matrix M(nobs, ncoeff);
   double elt;
   int p;
   typename EXPT_SET::const_iterator cit;
   typename std::vector< VAR_TUPLE >::const_iterator oit;
   int i,j;
   for (cit=expset.begin(), j=1; cit != expset.end() ; ++cit, ++j) //loop on exponent tuples
   {
      for(oit=obs_input.begin(), i=1; oit!=obs_input.end();++oit, ++i) //loop on observations
      {
         elt=1.0;
         for(int d=0;d<DIM;++d)
         {
            p = (*cit)[d];
            if (p != 0)
            {
               elt *= std::pow( (*oit)[d], p );
            }
         }
         M(i,j) = elt;
      }
   }

   NEWMAT::ColumnVector v(nobs);
   for(int o=0;o<nobs;++o)
   {
      v(o+1) = obs_output[o];
   }

   //perform solution
   NEWMAT::ColumnVector bfit = (M.t()*M).i()*M.t()*v;

   //compute RMS (optional, if rms non null)
   if (prms!=NULL)
   {
      NEWMAT::ColumnVector delta = M*bfit - v;
      *prms = std::sqrt( delta.SumSquare() / nobs);
   }

   //init polynom
   for (cit=expset.begin(), j=1; cit != expset.end() ; ++cit, ++j) //loop on exponent tuples
   {
      setMonom(*cit, bfit(j));
   }
   return true;
}

Member Data Documentation

theEpsilon

template<class T, int DIM = 1>

T ossimPolynom< T, DIM >::theEpsilon

protected

associate a scalar to each tuple of orders : monom

Definition at line 783 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::getEpsilon(), ossimPolynom< ossim_float64, 3 >::import(), ossimPolynom< ossim_float64, 3 >::isNegligible(), and ossimPolynom< ossim_float64, 3 >::operator=().

theMonoms

template<class T, int DIM = 1>

MONOM_MAP ossimPolynom< T, DIM >::theMonoms

protected

protected data members

Definition at line 781 of file ossimPolynom.h.

Referenced by ossimPolynom< ossim_float64, 3 >::addMonomQuick(), ossimPolynom< ossim_float64, 3 >::delMonom(), ossimPolynom< ossim_float64, 3 >::discardNullMonoms(), ossimPolynom< ossim_float64, 3 >::getCoeff(), ossimPolynom< ossim_float64, 3 >::getMonoms(), ossimPolynom< ossim_float64, 3 >::import(), ossimPolynom< ossim_float64, 3 >::isNull(), ossimPolynom< ossim_float64, 3 >::nullify(), ossimPolynom< ossim_float64, 3 >::operator=(), and ossimPolynom< ossim_float64, 3 >::setMonom().

---

The documentation for this class was generated from the following file:

• ossimPolynom.h