dfunc_evaluator.h

Go to the documentation of this file.

// Forward Differentiation Function Evaluator implementation -*- C++ -*-

// $Id: dfunc_evaluator.h 347 2008-02-06 15:20:01Z schilly $
// Copyright (C) 2001-2004 Hermann Schichl
// Copyright (C) 2005 Hermann Schichl, Ferenc Domes
//
// This file is part of the COCONUT API.  This library
// is free software; you can redistribute it and/or modify it under the
// terms of the Library GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library 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
// Library GNU General Public License for more details.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the Library GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the Library GNU General Public License.

#ifndef _DFUNC_EVALUATOR_H_
#define _DFUNC_EVALUATOR_H_

#include <coconut_config.h>
#include <evaluator.h>
#include <expression.h>
#include <model.h>
#include <eval_main.h>
#include <linalg.h>
#include <math.h>
#include <api_exception.h>

using namespace vgtl;

namespace coco {

typedef double (*func_evaluator)(const std::vector<double>* __x,
                                 const variable_indicator& __v);

template <typename _T, int DN>
struct dfunc_eval_rettype
{
  _T f[DN+1];
};

template <typename _T, int DN>
struct dfunc_eval_type
{
  unsigned int _0chnbase;
  _T x;
  std::vector<dfunc_eval_rettype<_T, DN> >* cache;
  const model* mod;
  void *p;
  _T d;
  dfunc_eval_rettype<_T, DN> r;
  unsigned int n;
};

template <typename _T, int DN>
class dfunc_eval :
  public cached_forward_evaluator_base<dfunc_eval_type<_T, DN>,expression_node,
                                       dfunc_eval_rettype<_T, DN>,
                                       expression_const_walker>
{
private:
  typedef dfunc_eval_rettype<_T, DN> retval;
  typedef cached_forward_evaluator_base<dfunc_eval_type<_T, DN>,expression_node,
                                        retval, expression_const_walker> _Base;
  typedef typename _Base::node_data_type node_data_type;

protected:
  bool is_cached(const node_data_type& __data)
  {
    if(__data.operator_type == EXPRINFO_LIN ||
       __data.operator_type == EXPRINFO_QUAD)
      return true;
    if(eval_data.cache && __data.n_parents > 1 && __data.n_children > 0 &&
       v_ind->match(__data.var_indicator()))
      return true;
    else
      return false;
  }

private:
  retval __power(double __coeff, retval __x, int __exp)
  {
    retval x;
    if(__exp == 0)
    {
      x.f[0] = _T(1.);
      for(int i=0; i <= DN; ++i)
        x.f[i] = _T(0.);
    }
    else
    {
      double k = __coeff*__x;
      switch(__exp)
      {
        case 1:
          return k;
          break;
        case 2:
          return k*k;
          break;
        case -1:
          return 1./k;
          break;
        case -2:
          return 1./(k*k);
          break;
        default:
          if(__exp & 1) // is odd
          {
            if(k < 0)
              return -std::pow(-k, __exp);
            else
              return std::pow(k, __exp);
          }
          else // even
            return std::pow(fabs(k), __exp);
          break;
      }
    }
  }

public:
  dfunc_eval(const _T& __x, const variable_indicator& __v,
            const model& __m, std::vector<retval>* __c) : _Base()
  {
    eval_data.x = __x;
    eval_data.cache = __c;
    eval_data.mod = &__m;
    eval_data.d = 0.;
    v_ind = &__v;
    eval_data.n = 0;
  }

  dfunc_eval(const dfunc_eval& __v) : _Base(__v) {}
  
  ~dfunc_eval() {}

  expression_const_walker short_cut_to(const expression_node& __data)
  { return eval_data.mod->node(0); }    // return anything
  

  void new_point(const _T& __x, const variable_indicator& __v)
  {
    eval_data.x = __x;
    v_ind = &__v;
  }

  void initialize() { return; }
  
  int initialize(const expression_node& __data)
  {
    int ret = 1;
    eval_data.n = 0;
#if 0
    if(__data.ev != NULL && (*__data.ev)[DFUNC_EVALUATOR] != NULL)
    // is there a short-cut evaluator defined?
    {
      eval_data.r = (*(func_evaluator)(*__data.ev)[DFUNC_EVALUATOR])(eval_data.x,
                                                                 *v_ind);
      return 0;     // don't perform the remaining graph walk
    }
#endif
    if(__data.node_num == eval_data._0chnbase)
    {
      eval_data.r.f[0] = eval_data.x;
      ret = 0;
    }
    else
    {
      switch(__data.operator_type)
      {
        case EXPRINFO_VARIABLE:
          // this must never occur!
          break;
        case EXPRINFO_SUM:
        case EXPRINFO_PROD:
        case EXPRINFO_MAX:
        case EXPRINFO_MIN:
        case EXPRINFO_INVERT:
        case EXPRINFO_DIV:
        case EXPRINFO_CONSTANT:
          eval_data.r.f[0] = __data.params.nd();
          break;
        case EXPRINFO_MONOME:
        case EXPRINFO_IN:
        case EXPRINFO_AND:
        case EXPRINFO_NOGOOD:
          eval_data.r.f[0] = 1.;
          break;
        case EXPRINFO_HISTOGRAM:
          eval_data.p = 
            (void*)new std::vector<unsigned int>(0,
                                           (__data.params.i().size()+1)>>1);
          eval_data.r.f[0] = 1.;
          break;
        case EXPRINFO_ALLDIFF:
          eval_data.p = (void*) new std::vector<double>;
          ((std::vector<double>*)eval_data.p)->reserve(__data.n_children);
          // no break!
        case EXPRINFO_MEAN:
        case EXPRINFO_IF:
        case EXPRINFO_OR:
        case EXPRINFO_NOT:
        case EXPRINFO_COUNT:
        case EXPRINFO_SCPROD:
        case EXPRINFO_NORM:
        case EXPRINFO_LEVEL:
          eval_data.r.f[0] = 0.;
          break;
        case EXPRINFO_DET:
        case EXPRINFO_PSD:
          // here construct square double matrix
          break;
        case EXPRINFO_COND:
        case EXPRINFO_FEM:
        case EXPRINFO_MPROD:
          // here construct rectangular double matrix
          break;
      }
      for(int i=1; i <= DN; ++i)
        eval_data.r.f[i] = _T(0.);
    }
    return ret;
  }

  void calculate(const expression_node& __data)
  {
#if 0
    if(__data.operator_type > 0)
    {
      eval_data.r = __data.f_evaluate(-1, __data.params.nn(), *eval_data.x,
                                      *v_ind, eval_data.r, 0, NULL);
    }
#endif
  }

  void retrieve_from_cache(const expression_node& __data)
  {
    eval_data.r = (*eval_data.cache)[__data.node_num];
  }

  int update(const retval& __rval)
  {
    eval_data.r = __rval;
    return 0;
  }

  int update(const expression_node& __data, const retval& __rval)
  {
    double __x;
    int ret = 0;
    if(__data.operator_type < 0)
    {
      switch(__data.operator_type)
      {
        case EXPRINFO_CONSTANT:
        case EXPRINFO_VARIABLE:
          break;
        case EXPRINFO_SUM:
        case EXPRINFO_MEAN:
          for(int i=0; i <= DN; ++i)
            eval_data.r.f[i] += __data.coeffs[eval_data.n]*__rval.f[i];
          eval_data.n++;
          break;
        case EXPRINFO_PROD:
          // the coefficients MUST be collected in __data.params.nd()
          {
            int _c[DN+1];
            int i, j, k, _co, _cn;
            memset(_c+1, 0, sizeof(int)*(DN+1));
            retval f = eval_data.r;
            for(i=0; i <= DN; _c[++i]=1, _co = 1)
            {
              eval_data.r.f[i] *= __rval.f[0];
              for(j=1, k=i-1; j <= i; ++j, --k, _co = _cn)
              {
                _cn = _c[j];
                eval_data.r.f[i] += _cn*__rval.f[j]*f[k];
                _c[j]+=_co;
              }
            }
          }
          break;
        case EXPRINFO_MONOME:
          eval_data.r *= __power(__data.coeffs[eval_data.n],
                        __rval, __data.params.n()[eval_data.n]);
          if(eval_data.r.f[0] == _T(0.)) ret=-1;
          eval_data.n++;
          break;
        case EXPRINFO_MAX:
          __x = __rval * __data.coeffs[eval_data.n++];
          if(__x > eval_data.r)
            eval_data.r = __x;
          break;
        case EXPRINFO_MIN:
          __x = __rval * __data.coeffs[eval_data.n++];
          if(__x < eval_data.r)
            eval_data.r = __x;
          break;
        case EXPRINFO_SCPROD:
          { double h = __data.coeffs[eval_data.n]*__rval;
            // if the number of children is odd, the list is implicitly
            // lengthened by one null element (i.e. implicitly shortened)
            if(eval_data.n & 1) // y_n
              eval_data.r += eval_data.d*h;
            else // x_n
              eval_data.d = h;
          }
          eval_data.n++;
          break;
        case EXPRINFO_NORM:
          if(__data.params.nd() == COCO_INF)
          {
            __x = fabs(__rval*__data.coeffs[eval_data.n]);
            if(__x > eval_data.r)
              eval_data.r = __x;
          }
          else
          { __x = fabs(__data.coeffs[eval_data.n]*__rval);
            eval_data.r += std::pow(__x, __data.params.nd());
          }
          eval_data.n++;
          if(eval_data.n == __data.n_children && __data.params.nd() != COCO_INF)
            eval_data.r = std::pow(eval_data.r, 1./__data.params.nd());
          break;
        case EXPRINFO_INVERT:
          eval_data.r /= __rval;
          break;
        case EXPRINFO_DIV:
          // this evaluator requires, that the second coefficient
          // is put into the first.
          if(eval_data.n == 0)
          {
            eval_data.r *= __rval;
            if(eval_data.r == 0) ret=-1;
          }
          else
            eval_data.r /= __rval;
          ++eval_data.n;
          break;
        case EXPRINFO_SQUARE:
          __x = __data.coeffs[0]*__rval+__data.params.nd();
          eval_data.r = __x*__x;
          break;
        case EXPRINFO_INTPOWER:
          eval_data.r = __power(__data.coeffs[0], __rval, __data.params.nn());
          break;
        case EXPRINFO_SQROOT:
          eval_data.r = std::sqrt(__data.coeffs[0]*__rval+__data.params.nd());
          break;
        case EXPRINFO_ABS:
          eval_data.r = fabs(__data.coeffs[0]*__rval+__data.params.nd());
          break;
        case EXPRINFO_POW:
          __x = __rval * __data.coeffs[eval_data.n];
          if(eval_data.n == 0)
            eval_data.r = __x+__data.params.nd();
          else
            eval_data.r = std::pow(eval_data.r, __x);
          ++eval_data.n;
          break;
        case EXPRINFO_EXP:
          eval_data.r = std::exp(__rval*__data.coeffs[0]+__data.params.nd());
          break;
        case EXPRINFO_LOG:
          eval_data.r = std::log(__rval*__data.coeffs[0]+__data.params.nd());
          break;
        case EXPRINFO_SIN:
          eval_data.r = std::sin(__rval*__data.coeffs[0]+__data.params.nd());
          break;
        case EXPRINFO_COS:
          eval_data.r = std::cos(__rval*__data.coeffs[0]+__data.params.nd());
          break;
        case EXPRINFO_ATAN2:
          __x = __rval * __data.coeffs[eval_data.n];
          if(eval_data.n == 0)
            eval_data.r = __x;
          else
            eval_data.r = std::atan2(eval_data.r, __x);
          ++eval_data.n;
          break;
        case EXPRINFO_GAUSS:
          __x = (__data.coeffs[0]*__rval-__data.params.d()[0])/
                                __data.params.d()[1];
          eval_data.r = std::exp(-__x*__x);
          break;
        case EXPRINFO_POLY:
          if(!__data.params.empty())
          {
            __x = __data.coeffs[eval_data.n]*__rval;
            eval_data.r = __data.params.d().back();
            for(int i = (int)__data.params.d().size()-2; i >= 0; --i)
            {
              eval_data.r *= __x;
              eval_data.r += __data.params.d()[i];
            }
          }
          else
            eval_data.r = 0;
          break;
        case EXPRINFO_LIN:
        case EXPRINFO_QUAD:
          /* we will never come here! */
          break;
        case EXPRINFO_IN:
          {
            __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(eval_data.r != -1 && i.contains(__x))
            {
              if(eval_data.r == 1 && (__x == i.inf() || __x == i.sup()))
                eval_data.r = 0;
            }
            else
            {
              eval_data.r = -1;
              ret = -1;
            }
          }
          eval_data.n++;
          break;
        case EXPRINFO_IF:
          __x = __rval * __data.coeffs[eval_data.n];
          if(eval_data.n == 0) // this is the condition
          {
            const interval& i(__data.params.ni());
            if(!i.contains(__x)) // go to else part
              ret = 1;           // don't walk the next child
          }
          else // this is either the if or the else part
          {
            eval_data.r = __x;
            ret = -1;            // don't continue walking
          }
          eval_data.n += ret+1;
          break;
        case EXPRINFO_AND:
          { __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(eval_data.r == 1 && !i.contains(__x))
            {
              eval_data.r = 0;
              ret = -1;         // result cannot become more false
            }
          }
          eval_data.n++;
          break;
        case EXPRINFO_OR:
          { __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(eval_data.r == 0 && i.contains(__x))
            {
              eval_data.r = 1;
              ret = -1;         // result cannot become more true
            }
          }
          eval_data.n++;
          break;
        case EXPRINFO_NOT:
          { __x = __data.coeffs[0]*__rval;
            const interval& i(__data.params.ni());
            if(i.contains(__x))
              eval_data.r = 0;
            else
              eval_data.r = 1;
          }
          break;
        case EXPRINFO_IMPLIES:
          { const interval& i(__data.params.i()[eval_data.n]);
            __x = __rval * __data.coeffs[eval_data.n];
            if(eval_data.n == 0)
            {
              if(!i.contains(__x))
              {
                eval_data.r = 1;
                ret = -1;        // ex falso quodlibet
              }
            }
            else // we only get here if first clause is true
              eval_data.r = i.contains(__x) ? 1 : 0;
            ++eval_data.n;
          }
          break;
        case EXPRINFO_COUNT:
          { __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(i.contains(__x))
              eval_data.r += 1;
          }
          eval_data.n++;
          break;
        case EXPRINFO_ALLDIFF:
          { __x = __data.coeffs[eval_data.n]*__rval;
            for(std::vector<double>::const_iterator _b =
                        ((std::vector<double>*)eval_data.p)->begin();
                _b != ((std::vector<double>*)eval_data.p)->end(); ++_b)
            {
              if(fabs(__x-*_b) <= __data.params.nd())
              {
                eval_data.r = 0;
                ret = -1;
                break;
              }
            }
            if(ret != -1)
              ((std::vector<double>*) eval_data.p)->push_back(__x);
          }
          eval_data.n++;
          if(eval_data.n == __data.n_children || ret == -1)
            delete (std::vector<double>*) eval_data.p;
          break;
        case EXPRINFO_HISTOGRAM:
          { __x = __data.coeffs[eval_data.n]*__rval;
            unsigned int ni = (__data.params.i().size()+1)>>1;
            for(unsigned int i=0; i<ni; i++)
            {
              if(__data.params.i()[i<<1].contains(__x))
                (*(std::vector<unsigned int>*)eval_data.p)[i]++;
            }
            eval_data.n++;
            if(eval_data.n == __data.n_children || ret == -1)
            {
              for(unsigned int i=0; i<ni; i++)
              {
                if(!__data.params.i()[i<<1].contains(
                      (*(std::vector<unsigned int>*)eval_data.p)[i]+0.))
                {
                  eval_data.r = 0;
                  break;
                }
              }
              delete (std::vector<unsigned int>*) eval_data.p;
            }
          }
          break;
        case EXPRINFO_LEVEL:
          { int h = (int)eval_data.r;
            __x = __data.coeffs[eval_data.n]*__rval;
            interval _h;

            if(h != INT_MAX)
            {
              while(h < __data.params.im().nrows())
              { 
                _h = __data.params.im()[h][eval_data.n];
                if(_h.contains(__x))
                  break;
                h++;
              }
              if(h == __data.params.im().nrows())
              {
                ret = -1;
                eval_data.r = INT_MAX;
              }
              else
                eval_data.r = h;
            }
          }
          eval_data.n++;
          break;
        case EXPRINFO_NEIGHBOR:
          if(eval_data.n == 0)
            eval_data.r = __data.coeffs[0]*__rval;
          else
          {
            double h = eval_data.r;
            eval_data.r = 0;
            __x = __data.coeffs[1]*__rval;
            for(unsigned int i = 0; i < __data.params.n().size(); i+=2)
            {
              if(h == __data.params.n()[i] && __x == __data.params.n()[i+1])
              {
                eval_data.r = 1;
                break;
              }
            }
          }
          eval_data.n++;
          break;
        case EXPRINFO_NOGOOD:
          {
            __x = __data.coeffs[eval_data.n]*__rval;
            if(eval_data.r == 0 || __data.params.n()[eval_data.n] != __x)
            {
              eval_data.r = 0;
              ret = -1;
            }
          }
          eval_data.n++;
          break;
        case EXPRINFO_EXPECTATION:
          throw nyi_exception("dfunc_evaluator: EXPECTATION");
          break;
        case EXPRINFO_INTEGRAL:
          throw nyi_exception("dfunc_evaluator: INTEGRAL");
          break;
        case EXPRINFO_LOOKUP:
        case EXPRINFO_PWLIN:
        case EXPRINFO_SPLINE:
        case EXPRINFO_PWCONSTLC:
        case EXPRINFO_PWCONSTRC:
          throw nyi_exception("dfunc_evaluator: Table Operations");
        case EXPRINFO_DET:
        case EXPRINFO_COND:
        case EXPRINFO_PSD:
        case EXPRINFO_MPROD:
        case EXPRINFO_FEM:
          throw nyi_exception("dfunc_evaluator: Matrix Operations");
          break;
        case EXPRINFO_RE:
        case EXPRINFO_IM:
        case EXPRINFO_ARG:
        case EXPRINFO_CPLXCONJ:
          throw nyi_exception("dfunc_evaluator: Complex Operations");
          break;
        case EXPRINFO_CMPROD:
        case EXPRINFO_CGFEM:
          throw nyi_exception("dfunc_evaluator: Const Matrix Operations");
        default:
          throw api_exception(apiee_evaluator,
              std::string("dfunc_evaluator: unknown function type ")+
              convert_to_str(__data.operator_type));
          break;
      }
    }
    else if(__data.operator_type > 0)
      // update the function arguments
      eval_data.r = __data.f_evaluate(eval_data.n++, __data.params.nn(),
                                      *eval_data.x, *v_ind, eval_data.r,
                                      __rval, NULL);
    // use the cache only if it is worthwhile
    if(eval_data.cache && __data.n_parents > 1 && __data.n_children > 0)
      (*eval_data.cache)[__data.node_num] = eval_data.r;
    return ret;
  }

  retval calculate_value(bool eval_all)
  {
    return eval_data.r;
  }
};

} // namespace coco

#endif /* _DFUNC_EVALUATOR_H_ */

---