OMKNumericPolator.h

Go to the documentation of this file.

/*
 * This file is part of openMask © INRIA, CNRS, Universite de Rennes 1 1993-2002, thereinafter the Software
 * 
 * The Software has been developped within the Siames Project. 
 * INRIA, the University of Rennes 1 and CNRS jointly hold intellectual property rights
 * 
 * The Software has been registered with the Agence pour la Protection des
 * Programmes (APP) under registration number IDDN.FR.001.510008.00.S.P.2001.000.41200
 *
 * This file may be distributed under the terms of the Q Public License
 * version 1.0 as defined by Trolltech AS of Norway and appearing in the file
 * LICENSE.QPL included in the packaging of this file.
 *
 * Licensees holding valid specific licenses issued by INRIA, CNRS or Université de Rennes 1 
 * for the software may use this file in accordance with that specific license 
 *
 */

#ifndef OMKPolatorNUMERIQUEHEADER
#define OMKPolatorNUMERIQUEHEADER

//template <class Type> class AbstractFifo;
#include <OMKPolator.h>
#include "OMKUserException.h"
#include "OMKTracer.h"
//------------------------------------------------------------------------
namespace OMK 
{
namespace Type
{

template <typename Type>
class NumericPolatorT : public Polator< SimpleTypeT< Type > >
{
 public:
  using Polator< SimpleTypeT<Type> >::get ;
  using Polator< SimpleTypeT<Type> >::getDate ;
  using Polator< SimpleTypeT<Type> >::getNumberOfPresentValues ;

 public:

  NumericPolatorT(void)  ;
  virtual ~NumericPolatorT(void) ;

  virtual const SimpleTypeT< Type > & 
    interpolate( SimpleTypeT< Type > & resultPlaceHolder,
                const int interprecisionLevel,
                const Date & dateNeeded,
                const Date & dateAfter,
                const SimpleTypeT<Type> & valueAfter,
                const Date & dateBefore,
                int offsetToMostRecentOfDateBefore ) const ;
  
  virtual const SimpleTypeT< Type > & 
    extrapolate( SimpleTypeT< Type > & resultPlaceHolder,
                const int requestedPrecisionLevel,
                const Date & t,
                const Date & tIndice ) const ;
  virtual const SimpleTypeT< Type > & 
    antepolate( SimpleTypeT< Type > & resultPlaceHolder,
                const int requestedPrecisionLevel,
                const Date & t,
                const Date & tIndice,
                unsigned int indice ) const ;
 private:
  static Type convert( double v ) ;
 public:
  static Type linearInterpolate( long t,
                                 long t1, Type v1,
                                 long t0, Type v0 ) 
  {
    return linearExtrapolate( t, t1,  v1, t0, v0 ) ;
  }

  static Type quadraticInterpolate( long t,
                                    long t2, Type v2,
                                    long t1, Type v1,
                                    long t0, Type v0 ) 
  {
    return quadraticExtrapolate( t, t2, v2, t1, v1, t0, v0 ) ;
  }


  static Type cubicInterpolate( long t,
                                long t3, Type v3,
                                long t2, Type v2,
                                long t1, Type v1,
                                long t0, Type v0 ) 
  {
    return cubicExtrapolate( t, t3, v3, t2, v2, t1, v1, t0, v0 ) ;
  }



  static Type linearExtrapolate( long t,
                                 long t1, Type v1,
                                 long t0, Type v0 )
  {
    double d0 = t1 - t0 ; OMASSERTM ( d0, "Will divide by 0" ) ;
    double f0 = t - t0 ; 
    // formula :  ((t - t0) / d0) * v1  +((t1 - t) / d0)* v0 ;
    return convert( ( v1 - v0 ) * ( f0 / d0 ) + v0 ) ;
  }


  static Type quadraticExtrapolate( long t,
                                    long t2, Type v2,
                                    long t1, Type v1,
                                    long t0, Type v0 ) 
  {
    double d0 = t1 - t0 ;  OMASSERTM( d0, "Will divide by 0" ) ;
    double d1 = t2 - t0 ;  OMASSERTM( d1, "Will divide by 0" ) ;
    double d2 = t1 - t2 ;  OMASSERTM( d2, "Will divide by 0" ) ;
    
    double k0 = ( t - t1 ) / d1 ;
    double k1 = ( t - t2 ) / d0 ;
    double k2 = ( t - t0 ) / d2 ;
    // formula : (v1*(k2*k1)) - (v2*(k2*k0))  + ((k0*k1)*v0)
    return convert( ( ( v1 - v0 ) * k1 - ( v2 - v0 ) * k0 ) * k2 + v0 ) ; 
  }


  static Type cubicExtrapolate( long t,
                                long t3, Type v3,
                                long t2, Type v2,
                                long t1, Type v1,
                                long t0, Type v0 ) 
  {
    double d0 = t0 - t1 ;  OMASSERTM( d0, "Will divide by 0" ) ;
    double d1 = t2 - t0 ;  OMASSERTM( d1, "Will divide by 0" ) ;
    double d2 = t0 - t3 ;  OMASSERTM( d2, "Will divide by 0" ) ;
    double d3 = t1 - t2 ;  OMASSERTM( d3, "Will divide by 0" ) ;
    double d4 = t1 - t3 ;  OMASSERTM( d4, "Will divide by 0" ) ;
    double d5 = t3 - t2 ;  OMASSERTM( d5, "Will divide by 0" ) ;
    
    double f0 = t - t0 ;
    double f1 = t - t1 ;
    double f2 = t - t2 ;
    double f3 = t - t3 ;

    double k0 = f1 * ( f2 / d2 ) ;
    double k1 = f0 * ( f3 / d3 ) ;
    //  formula : v3*((f0/d5)*(k0/d4))+v2*((f1/d1)*(k1/d5))-v1*((k1/d0)*(f2/d4))-v0*((k0/d0)*(f3/d1))
    return convert( ( v3 - v0 ) * ( ( f0 / d5 ) * ( k0 / d4 ) )
                 + ( v2 - v0 ) * ( ( f1 / d1 ) * ( k1 / d5 ) )
                 - ( v1 - v0 ) * ( ( k1 / d0 ) * ( f2 / d4 ) )
                 + v0 ) ; 
    }
};
template <>
inline int NumericPolatorT<int>::convert( double v ) 
{
  return (int)v ;
}
template <>
inline float NumericPolatorT<float>::convert( double v ) 
{
  return (float)v ;
}
template <>
inline double NumericPolatorT<double>::convert( double v ) 
{
  return v ;
}
/*
template < typename Type >
inline Type NumericPolatorT< Type >::convert( double v ) 
{
  return (Type)v ;
}
*/

//---------------------------------------------------------------
template <typename Type>
inline NumericPolatorT<Type>::NumericPolatorT() : Polator< SimpleTypeT<Type> > (3)
{
}
  
//---------------------------------------------------------------
template <typename Type>
inline NumericPolatorT<Type>::~NumericPolatorT()
{
}
//------------------------------------------------------------------------

template <typename Type>
inline const SimpleTypeT<Type> & 
NumericPolatorT<Type>::interpolate( SimpleTypeT< Type >& resultPlaceHolder,
                                       const int interprecisionLevel,
                                       const Date & dateNeeded,
                                       const Date & dateAfter,
                                       const SimpleTypeT<Type> & valueAfter,
                                       const Date & dateBefore,
                                       int offsetToMostRecentOfDateBefore ) const 
{
   
#ifdef _DEBUGTYPEUTIL
  cout << "*********************************************" << endl ;
  cout << "NumericPolatorT<Type>::interpolate" << endl;
  cout << "interpolate level : " << interprecisionLevel <<endl ;
  cout << "Date dateAfter  <----- : " << dateAfter << endl ;
  cout << "Date dateNeeded <----- : " << dateNeeded << endl ;
  cout << "Date dateBefore <----- : " << dateBefore << endl ;
  cout << "Type : valueAfter " << valueAfter << endl ;
  cout << "int : offsetToMostRecentOfDateBefore " << offsetToMostRecentOfDateBefore << endl ;
  cout << "*********************************************" << endl ;
#endif
    
  switch ( interprecisionLevel ) 
  {
  case PolatorNT::Constant : 
  case PolatorNT::ConstantContinuous :
    resultPlaceHolder = ( ( dateAfter - dateNeeded ) <= ( dateNeeded - dateBefore ) ) ?
                          valueAfter : get (offsetToMostRecentOfDateBefore) ;
    break ;
  case PolatorNT::Linear : 
  case PolatorNT::LinearContinuous :
    resultPlaceHolder = 
      linearInterpolate( dateNeeded, 
                         dateAfter, (Type)valueAfter, 
                         dateBefore, (Type)get( offsetToMostRecentOfDateBefore ) ) ;
    break ;
  case PolatorNT::Quadratic :
  case PolatorNT::QuadraticContinuous :
    resultPlaceHolder = 
      quadraticInterpolate( dateNeeded, 
                            dateAfter, (Type)valueAfter, 
                            dateBefore, (Type)get (offsetToMostRecentOfDateBefore) ,
                            getDate( offsetToMostRecentOfDateBefore + 1 ), 
                            (Type)get( offsetToMostRecentOfDateBefore + 1 ) ) ;
    break ;
  case PolatorNT::Cubic :
  case PolatorNT::CubicContinuous :
    resultPlaceHolder = 
      cubicInterpolate( dateNeeded, 
                        dateAfter, (Type)valueAfter, 
                        dateBefore, (Type)get( offsetToMostRecentOfDateBefore ),
                        getDate( offsetToMostRecentOfDateBefore + 1 ), 
                        (Type)get( offsetToMostRecentOfDateBefore + 1 ), 
                        getDate( offsetToMostRecentOfDateBefore + 2 ), 
                        (Type)get( offsetToMostRecentOfDateBefore + 2 ) ) ;
    break ;
  default :
    OMTRACEID( OMK_DEBUG_OMK_TYPE, "NumericPolatorT<Type>::interpolate : unknown polation level, using PolatorNT::Constant" ) ;
    resultPlaceHolder = ( ( dateAfter - dateNeeded ) <= ( dateNeeded - dateBefore ) ) ?
                          valueAfter : get( offsetToMostRecentOfDateBefore ) ;
    break ;
  }
  return resultPlaceHolder ;
}

//------------------------------------------------------------------------
template <typename Type>
inline 
const SimpleTypeT<Type>& 
NumericPolatorT<Type>::extrapolate( SimpleTypeT<Type> & resultPlaceHolder,
                                       const int requestedPrecisionLevel,
                                       const Date & t,
                                       const Date & tIndice )  const 
{
  //extrapolate à la sauce David Margery
  //il s'agit d'une extrapolation par morceaux
  //on cherche à trouver la valeur à t
#ifdef _DEBUGPOLATION
  cerr << "NumericPolatorT<Type>::extrapolate (" <<requestedPrecisionLevel<<", "<<t<<", "<<tIndice<<")"<<endl;
#endif
  switch (requestedPrecisionLevel) 
  {
  case PolatorNT::Constant : 
    resultPlaceHolder = get( 0 ) ;
    break ;
  case PolatorNT::ConstantContinuous : 
    resultPlaceHolder = get( 0 ) ;
    if( 1 < getNumberOfPresentValues() )
    {
      resultPlaceHolder.setValue( 
        linearInterpolate( t, tIndice, (Type)get( 1 ), 
                           2 * tIndice - getDate( 1 ), (Type)get( 0 ) ) ) ;
    }
    break ;
  case PolatorNT::Linear : 
    resultPlaceHolder.setValue( 
      linearExtrapolate( t, tIndice, (Type)get( 0 ), 
                         getDate( 1 ), (Type)get( 1 ) ) ) ;
    break ; 
  case PolatorNT::LinearContinuous : 
    resultPlaceHolder.setValue( 
      linearExtrapolate( t, tIndice, (Type)get( 0 ), 
                         getDate( 1 ), (Type)get( 1 ) ) ) ;
    if( 2 < getNumberOfPresentValues() )
    {
      Type oldResult =
        linearExtrapolate( t, getDate( 1 ), (Type)get( 1 ), 
                           getDate( 2 ), (Type)get( 2 ) ) ;
      resultPlaceHolder.setValue( 
        linearInterpolate( t, tIndice, oldResult, 
                           2 * tIndice - getDate( 1 ), (Type)resultPlaceHolder ) ) ;
    }
    break ; 

  case PolatorNT::Quadratic :
    resultPlaceHolder.setValue( 
      quadraticExtrapolate( t, tIndice, (Type)get( 0 ), 
                            getDate( 1 ), (Type)get( 1 ), 
                            getDate( 2 ), (Type)get( 2 ) ) ) ;
    break;
  case PolatorNT::QuadraticContinuous :
    resultPlaceHolder.setValue( 
      quadraticExtrapolate( t, tIndice, (Type)get( 0 ), 
                            getDate( 1 ), (Type)get( 1 ), 
                            getDate( 2 ), (Type)get( 2 ) ) ) ;

    if( getNumberOfPresentValues() > 3 )
    {
      Type oldResult= 
        quadraticExtrapolate( t, getDate( 1 ), (Type)get( 1 ), 
                              getDate( 2 ), (Type)get( 2 ), 
                              getDate( 3 ), (Type)get( 3 ) ) ;
      resultPlaceHolder.setValue( 
        linearInterpolate( t, tIndice, oldResult, 
                           2 * tIndice - getDate( 1 ), (Type)resultPlaceHolder ) ) ;
    }
    break;

  case PolatorNT::Cubic :
    resultPlaceHolder.setValue( cubicExtrapolate( t, tIndice, (Type)get( 0 ), 
                                          getDate( 1 ), (Type)get( 1 ), 
                                          getDate( 2 ), (Type)get( 2 ), 
                                          getDate( 3 ), (Type)get( 3 ) ) ) ;
    break;

  case PolatorNT::CubicContinuous :
    resultPlaceHolder.setValue( cubicExtrapolate( t, tIndice, (Type)get( 0 ),
                                          getDate( 1 ), (Type)get( 1 ),
                                          getDate( 2 ), (Type)get( 2 ),
                                          getDate( 3 ), (Type)get( 3 ) ) ) ;
    if( 4 < getNumberOfPresentValues() )
      {
        Type oldResult= 
          cubicExtrapolate( t, getDate( 1 ), (Type)get( 1 ),
                            getDate( 2 ), (Type)get( 2 ),
                            getDate( 3 ), (Type)get( 3 ),
                            getDate( 4 ), (Type)get( 4 ) ) ;
        resultPlaceHolder.setValue( 
          linearInterpolate( t, tIndice, oldResult, 
                             2*tIndice - getDate( 1 ), (Type)resultPlaceHolder ) ) ;
      }

    break;
  default :
    OMTRACEID( OMK_DEBUG_OMK_TYPE, "NumericPolatorT<Type>::interpolate : unknown polation level, using PolatorNT::Constant" ) ;
    resultPlaceHolder = get( 0 ) ;
  }
  return resultPlaceHolder ;
}
//------------------------------------------------------------------------
template <typename Type>
inline 
const SimpleTypeT<Type> & 
NumericPolatorT<Type>::antepolate ( SimpleTypeT< Type >& resultPlaceHolder,
                                       const int requestedPrecisionLevel,
                                       const Date& t,
                                       const Date& tIndice, 
                                       unsigned int indice ) const 
{
   
#ifdef _DEBUGTYPEUTIL
  cout << "*********************************************" << endl ;
  cout << "NumericPolatorT<Type>::antepolate" << endl;
  cout << "Date t       <----- : " << t << endl ;
  cout << "Date tIndice <----- : " << tIndice << endl ;
  cout << "polation level : " << requestedPrecisionLevel << endl ;
  cout << "index " << indice << endl ;
  cout << "*********************************************" << endl ;
#endif

  switch (requestedPrecisionLevel) 
  {
  case PolatorNT::Constant : 
  case PolatorNT::ConstantContinuous : 
    resultPlaceHolder = get( indice ) ;
    break ;
  case PolatorNT::Linear : 
  case PolatorNT::LinearContinuous : 
    resultPlaceHolder = 
      linearExtrapolate( t, tIndice, (Type)get( indice ), 
                         getDate( indice - 1 ), (Type)get( indice - 1 ) ) ;
    break ; 
  case PolatorNT::Quadratic :
  case PolatorNT::QuadraticContinuous :
    resultPlaceHolder = 
      quadraticExtrapolate( t, tIndice, (Type)get( indice ), 
                            getDate( indice - 1 ), (Type)get( indice - 1 ), 
                            getDate( indice - 2 ), (Type)get( indice - 2 ) ) ;
    break;
  case PolatorNT::Cubic :
  case PolatorNT::CubicContinuous :
    resultPlaceHolder = 
      cubicExtrapolate( t, tIndice, (Type)get( indice ), 
                        getDate( indice - 1 ), (Type)get( indice - 1 ), 
                        getDate( indice - 2 ), (Type)get( indice - 2 ), 
                        getDate( indice - 3 ), (Type)get( indice - 3 ) ) ;
    break;
  default :
    OMTRACEID( OMK_DEBUG_OMK_TYPE, "NumericPolatorT<Type>::interpolate : unknown polation level, using PolatorNT::Constant" ) ;
    resultPlaceHolder = get(indice) ;
  }
  return resultPlaceHolder ;
}
} // namespace Type
} // namespace OMK


#endif

---