/// @addtogroup dynsys
/// @{

/////////////////////////////////////////////////////////////////////////////
/// @file CnTaylor_remainder.hpp
///
/// @author Daniel Wilczak
/////////////////////////////////////////////////////////////////////////////

// Copyright (C) 2000-2012 by the CAPD Group.
//
// This file constitutes a part of the CAPD library,
// distributed under the terms of the GNU General Public License.
// Consult  http://capd.ii.uj.edu.pl/ for details.

#ifndef _CAPD_DYNSYS_CNTAYLOR_REMAINDER_HPP_
#define _CAPD_DYNSYS_CNTAYLOR_REMAINDER_HPP_

namespace capd{
namespace dynsys{

template<typename MapType, typename StepControlType,typename CurveType>
typename CnTaylor<MapType,StepControlType,CurveType>::VectorType
CnTaylor<MapType,StepControlType,CurveType>::Remainder(const ScalarType& t, const VectorType &iv, VectorType& o_enc)
{
  o_enc = this->enclosure(t,iv);
  const static ScalarType I(TypeTraits<ScalarType>::zero().leftBound(),TypeTraits<ScalarType>::one().rightBound());
  this->computeRemainderCoefficients(t + I*this->m_step,o_enc);
  size_type r=this->getOrder()+1;

  ScalarType fac = power(this->m_step,r);
  VectorType result(this->dimension(),true);
  for(size_type i=0;i<this->dimension();++i)
    result[i] = this->remainderCoefficient(i,r)*fac;

  return result;
}

//###########################################################//

template<typename MapType, typename StepControlType,typename CurveType>
void CnTaylor<MapType,StepControlType,CurveType>::JacRemainder(
      const ScalarType& t,
      const VectorType& vecEnclosure,
      const MatrixType& jacEnclosure,
      VectorType& Remainder,
      MatrixType& jacRemainder
    )
{
  const static ScalarType I(TypeTraits<ScalarType>::zero().leftBound(),TypeTraits<ScalarType>::one().rightBound());
  size_type r=this->getOrder()+1;
  this->computeRemainderCoefficients(t + I*this->m_step,vecEnclosure,jacEnclosure);
  ScalarType fac = power(this->m_step,r);

  for(size_type i=0;i<this->dimension();++i)
  {
    Remainder[i] = fac*this->remainderCoefficient(i,r);
    for(size_type j=0;j<this->dimension();++j)
      jacRemainder(i+1,j+1) = fac*this->remainderCoefficient(i,j,r);
  }
}

// ####################################################################

template<typename MapType,typename StepControlT, typename CurveT>
void CnTaylor<MapType,StepControlT,CurveT>::c2Remainder(
      const VectorType& Enc,
      const MatrixType& jacEnc,
      const HessianType& hessianEnc,
      VectorType& o_Rem,
      MatrixType& o_jacRem,
      HessianType& o_hessianRem
  )
{
  size_type i,j,c;
  // set initial condition
  for(i=0;i<this->dimension();++i)
  {
    this->remainderCoefficient(i,0) = Enc[i];
    for(j=0;j<this->dimension();++j)
    {
      this->remainderCoefficient(i,j,0) = jacEnc(i+1,j+1);
      for(c=j;c<this->dimension();++c)
        this->remainderCoefficient(i,j,c,0) = hessianEnc(i,j,c);
    }
  }

  this->m_vField->computeODECoefficients(this->getRemainderCoefficients(),2,this->getOrder()+1);

  ScalarType fac = power(this->m_step,this->getOrder()+1);
  for(i=0;i<this->dimension();++i)
  {
    o_Rem[i] = fac*this->remainderCoefficient(i,this->getOrder()+1);
    for(j=0;j<this->dimension();++j)
    {
      o_jacRem(i+1,j+1) = fac*this->remainderCoefficient(i,j,this->getOrder()+1);
      for(c=j;c<this->dimension();++c)
        o_hessianRem(i,j,c) = fac*this->remainderCoefficient(i,j,c,this->getOrder()+1);
    }
  }
}



}} // namespace capd::dynsys

#endif // _CAPD_DYNSYS_CNTAYLOR_REMAINDER_HPP_

/// @}
