---

der_evaluator.h

Go to the documentation of this file.

// Derivative Evaluator implementation -*- C++ -*-

// Copyright (C) 2001-2003 Hermann Schichl
//
// 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 _DER_EVALUATOR_H_
#define _DER_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>

using namespace vgtl;

typedef bool (*prep_d_evaluator)();
typedef double (*func_d_evaluator)(const std::vector<double>* __x,
                  const variable_indicator& __v,
                  std::vector<double>& __d_data);
typedef std::vector<double>& (*der_evaluator)(const std::vector<double>& __d_dat,
                                         const variable_indicator& __v);

class prep_d_eval : public
  cached_forward_evaluator_base<std::vector<std::vector<double> >*, expression_node, bool, model::walker>
{
private:
  typedef cached_forward_evaluator_base<std::vector<std::vector<double> >*,
                                    expression_node,bool,model::walker> _Base;

public:
  prep_d_eval(std::vector<std::vector<double> >& __d,
              unsigned int _num_of_nodes)
  {
    eval_data = &__d;
    if((*eval_data).size() < _num_of_nodes)
      (*eval_data).insert((*eval_data).end(),
                      _num_of_nodes-(*eval_data).size(), std::vector<double>());
  }

  prep_d_eval(const prep_d_eval& __x) { eval_data = __x.eval_data; }

  ~prep_d_eval() {}

  void initialize() { return; }

  bool is_cached(const expression_node& __data)
  {
    return (*eval_data)[__data.node_num].size() > 0;
  }

  void retrieve_from_cache(const expression_node& __data) { return; }

  int initialize(const expression_node& __data)
  {
    (*eval_data)[__data.node_num].insert((*eval_data)[__data.node_num].end(),
        __data.n_children, 0);
    return 1;
  }

  void calculate(const expression_node& __data) { return; }

  int update(bool __rval) { return 0; }

  int update(const expression_node& __data, bool __rval)
        { return 0; }

  bool calculate_value(bool eval_all) { return true; }
};

// function evaluation at double argument + preparation of derivative data
// for later backwards evaluation

struct func_d_eval_type
{
  const std::vector<double>* x;
  std::vector<double>* f_cache;
  std::vector<std::vector<double> >* d_data;
  const model* mod;
  union { void* p; double d; } u;
  double r;
  unsigned int n, 
           info; // needed for derivative of max, min
};

class func_d_eval : public
  cached_forward_evaluator_base<func_d_eval_type, expression_node,
                                double, model::walker>
{
private:
  typedef cached_forward_evaluator_base<func_d_eval_type,expression_node,
                                        double, model::walker> _Base;

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.f_cache && __data.n_parents > 1 && __data.n_children > 0 &&
       v_ind->match(__data.var_indicator()))
      return true;
    else
      return false;
  }

private:
  double __power(double __coeff, double __x, int __exp)
  {
    if(__exp == 0)
      return 1.;
    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 -pow(-k, __exp);
            else
              return pow(k, __exp);
          }
          else // even
            return pow(fabs(k), __exp);
          break;
      }
    }
  }

  void __calc_max(double h, const expression_node& __data)
  {
    if(h >= eval_data.r)
    {
      (*eval_data.d_data)[__data.node_num][eval_data.info] = 0;
      if(h > eval_data.r)
        (*eval_data.d_data)[__data.node_num][eval_data.n] =
                          __data.coeffs[eval_data.n];
      else
        // the derivative is here not defined, anyway
        (*eval_data.d_data)[__data.node_num][eval_data.n] = 0;
      eval_data.r = h;
      eval_data.info = eval_data.n;
    }
    else
    {
      (*eval_data.d_data)[__data.node_num][eval_data.n] = 0;
    }
  }

public:
  func_d_eval(const std::vector<double>& __x, const variable_indicator& __v,
              const model& __m, std::vector<std::vector<double> >& __d,
              std::vector<double>* __c) : _Base()
  {
    eval_data.x = &__x;
    eval_data.f_cache = __c;
    eval_data.mod = &__m;
    eval_data.d_data = &__d;
    eval_data.n = 0;
    eval_data.r = 0;
    eval_data.u.d = 0;
    v_ind = &__v;
  }

  func_d_eval(const func_d_eval& __x) : _Base(__x) {}

  ~func_d_eval() {}

  model::walker short_cut_to(const expression_node& __data)
      { return eval_data.mod->node(0); }  // return anything - not needed

  void new_point(const std::vector<double>& __x, const variable_indicator& __v)
  {
    eval_data.x = &__x;
    v_ind = &__v;
  }

  void initialize() { return; }

  int initialize(const expression_node& __data)
  {
    eval_data.n = 0;
    if(__data.ev != NULL && (*__data.ev)[FUNC_D_EVALUATOR] != NULL)
    // is there a short-cut evaluator defined?
    {
      eval_data.r =
          (*(func_d_evaluator)(*__data.ev)[FUNC_D_EVALUATOR])(eval_data.x,
                               *v_ind, (*eval_data.d_data)[__data.node_num]);
      return 0;     // don't perform the remaining graph walk
    }
    else
    {
      switch(__data.operator_type)
      {
        case EXPRINFO_MAX:
        case EXPRINFO_MIN:
          eval_data.info = 0;
          // no break == continue
        case EXPRINFO_SUM:
        case EXPRINFO_PROD:
        case EXPRINFO_INVERT:
          eval_data.r = __data.params.nd();
          break;
        case EXPRINFO_IN:
        case EXPRINFO_AND:
        case EXPRINFO_NOGOOD:
          eval_data.r = 1.;
          break;
        case EXPRINFO_ALLDIFF:
          eval_data.u.p = (void*) new std::vector<double>;
          ((std::vector<double>*)eval_data.u.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_LEVEL:
          eval_data.r = 0.;
          break;
        case EXPRINFO_NORM:
          eval_data.info = 0;
          eval_data.r = 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;
      }
      return 1;
    }
  }

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

  void retrieve_from_cache(const expression_node& __data)
  {
    // the parallel der_cache MUST be initialized, too
    if(__data.operator_type == EXPRINFO_LIN)
      eval_data.r = linalg_dot(eval_data.mod->lin[__data.params.nn()],
                               *eval_data.x,0.);
    else if(__data.operator_type == EXPRINFO_QUAD)
    {
      std::vector<double> irslt = *eval_data.x;
      unsigned int r = __data.params.m().nrows();

      // compute x^T A = irslt
      irslt.push_back(0);
      linalg_matvec(__data.params.m(), irslt, irslt);
      irslt.pop_back();
      // compute b^T x
      eval_data.r = linalg_dot(__data.params.m()[r-1], *eval_data.x, 0.);
      // compute x^T A x + b^T x
      eval_data.r += linalg_dot(irslt,*eval_data.x,0.);
      // compute x^T A + 0.5 b (half the derivative)
      linalg_add(linalg_scale(__data.params.m()[r-1], 0.5), irslt, 
                 (*eval_data.d_data)[__data.node_num]);
    }
    else
      eval_data.r = (*eval_data.f_cache)[__data.node_num];
  }

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

  int update(const expression_node& __data, const double& __rval)
  {
    int ret = 0;
    double __x;
    if(__data.operator_type < 0)
    {
      switch(__data.operator_type)
      {
        case EXPRINFO_CONSTANT:
          eval_data.r = __data.params.nd();
          // don't care about d_data, CONSTANTS are ALWAYS leafs
          break;
        case EXPRINFO_VARIABLE:
          eval_data.r = (*eval_data.x)[__data.params.nn()];
          // don't care about d_data, VARIABLES are ALWAYS leafs
          break;
        case EXPRINFO_SUM:
        case EXPRINFO_MEAN:
          { double h = __data.coeffs[eval_data.n];
            eval_data.r += h*__rval;
            (*eval_data.d_data)[__data.node_num][eval_data.n++] = h;
          }
          break;
        case EXPRINFO_PROD:
          if(eval_data.n == 0)
          {
            eval_data.r *= __rval;
            (*eval_data.d_data)[__data.node_num][0] = __data.params.nd();
          }
          else
          {
            (*eval_data.d_data)[__data.node_num][eval_data.n] = eval_data.r;
            eval_data.r *= __rval;
            for(int i = eval_data.n-1; i >= 0; i--)
              (*eval_data.d_data)[__data.node_num][i] *= __rval;
          }
          ++eval_data.n;
          break;
        case EXPRINFO_MONOME:
          if(eval_data.n == 0)
          {
            int n = __data.params.n()[0];
            if(n != 0)
            {
              __x = __power(__data.coeffs[0], __rval, n-1)*__data.coeffs[0];
              eval_data.r = __x*__rval;
              (*eval_data.d_data)[__data.node_num][0] = n*__x;
            }
            else
            {
              (*eval_data.d_data)[__data.node_num][0] = 0;
              eval_data.r = 1.;
            }
          }
          else
          {
            int n = __data.params.n()[eval_data.n];
            if(n != 0)
            {
              __x = __power(__data.coeffs[eval_data.n], __rval, n-1)*
                        __data.coeffs[eval_data.n];
              (*eval_data.d_data)[__data.node_num][eval_data.n] =
                                                           eval_data.r*n*__x;
              __x *= __rval;
              eval_data.r *= __x;
              for(int i = eval_data.n-1; i >= 0; i--)
                (*eval_data.d_data)[__data.node_num][i] *= __x;
            }
            else
              (*eval_data.d_data)[__data.node_num][eval_data.n] = 0;
          }
          ++eval_data.n;
          break;
        case EXPRINFO_MAX:
          __calc_max(__rval * __data.coeffs[eval_data.n], __data);
          ++eval_data.n;
          break;
        case EXPRINFO_MIN:
          { double h = __rval * __data.coeffs[eval_data.n];
            if(h <= eval_data.r)
            {
              (*eval_data.d_data)[__data.node_num][eval_data.info] = 0;
              if(h < eval_data.r)
                (*eval_data.d_data)[__data.node_num][eval_data.n] =
                                  __data.coeffs[eval_data.n];
              else
                // the derivative is here not defined, anyway
                (*eval_data.d_data)[__data.node_num][eval_data.n] = 0;
              eval_data.r = h;
              eval_data.info = eval_data.n;
            }
            else
            {
              (*eval_data.d_data)[__data.node_num][eval_data.n] = 0;
            }
          }
          ++eval_data.n;
          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.u.d*h;
              (*eval_data.d_data)[__data.node_num][eval_data.n] =
                                    eval_data.u.d*__data.coeffs[eval_data.n-1];
              (*eval_data.d_data)[__data.node_num][eval_data.n-1] =
                                                  h*__data.coeffs[eval_data.n];
            }
            else // x_n
              eval_data.u.d = h;
          }
          eval_data.n++;
          break;
        case EXPRINFO_NORM:
          if(__data.params.nd() == INFINITY)
            __calc_max(fabs(__rval * __data.coeffs[eval_data.n]), __data);
          else
          {
            double h = __data.coeffs[eval_data.n]*fabs(__rval);
            double O = pow(h, __data.params.nd()-1);
            eval_data.r += O*h;
            (*eval_data.d_data)[__data.node_num][eval_data.n] = O;
          }             
          eval_data.n++;
          if(eval_data.n == __data.n_children &&
             __data.params.nd() != INFINITY)
          {
            double h = pow(eval_data.r,1./(__data.params.nd())-1.);
            for(unsigned int i = 0; i < eval_data.n; ++i)
              (*eval_data.d_data)[__data.node_num][eval_data.n] *= h;
            eval_data.r = pow(eval_data.r,1./(__data.params.nd()));
          }
          break;
        case EXPRINFO_INVERT:
          { double h = 1/__rval;
            eval_data.r *= h;
            (*eval_data.d_data)[__data.node_num][0] = -__data.params.nd()*h*h;
          }
          break;
        case EXPRINFO_DIV:  // has two children
          if(eval_data.n++ == 0)
            eval_data.r = __rval;
          else
          {
            double h = 1/__rval;
            eval_data.r *=
                (*eval_data.d_data)[__data.node_num][0] = __data.params.nd()*h;
            (*eval_data.d_data)[__data.node_num][1] = -eval_data.r*h;
          }
          break;
        case EXPRINFO_SQUARE:
          { double h = __data.coeffs[0]*__rval+__data.params.nd();
            eval_data.r = h*h;
            (*eval_data.d_data)[__data.node_num][0] = 2*h*__data.coeffs[0];
          }
          break;
        case EXPRINFO_INTPOWER:
          { int hl = __data.params.nn();
            if(hl == 0)
            {
              eval_data.r = 1;
              (*eval_data.d_data)[__data.node_num][0] = 0;
            }
            else
            {
              double kl = __data.coeffs[0]*__rval;
              switch(hl)
              {
                case 1:
                  eval_data.r = kl;
                  (*eval_data.d_data)[__data.node_num][0] = __data.coeffs[0];
                  break;
                case 2:
                  eval_data.r = kl*kl;
                  (*eval_data.d_data)[__data.node_num][0] =
                        2*kl*__data.coeffs[0];
                  break;
                case -1:
                  { double h = 1/kl;
                    eval_data.r = h;
                    (*eval_data.d_data)[__data.node_num][0] =
                                -h*h*__data.coeffs[0];
                  }
                  break;
                case -2:
                  { double h = 1/kl;
                    double k = h*h;
                    eval_data.r = k;
                    (*eval_data.d_data)[__data.node_num][0] =
                                                -2*h*k*__data.coeffs[0];
                  }
                  break;
                default:
                  { double h;
                    if(hl & 1) // is odd
                      h = pow(fabs(kl), hl-1);
                    else // even
                    {
                      if(kl < 0)
                        h = -pow(-kl, hl-1);
                      else
                        h = pow(kl, hl-1);
                    }
                    eval_data.r = h*kl;
                    (*eval_data.d_data)[__data.node_num][0] =
                                        hl*h*__data.coeffs[0];
                  }
                  break;
              }
            }
          }
          break;
        case EXPRINFO_SQROOT:
          { double h = sqrt(__data.coeffs[0]*__rval+__data.params.nd());
            eval_data.r = h;
            (*eval_data.d_data)[__data.node_num][0] = 0.5*__data.coeffs[0]/h;
          }
          break;
        case EXPRINFO_ABS:
          { double h = __data.coeffs[0]*__rval+__data.params.nd();
            eval_data.r = fabs(h);
            (*eval_data.d_data)[__data.node_num][0] =
                    h > 0 ? __data.coeffs[0] : (h < 0 ? -__data.coeffs[0] : 0);
          }
          break;
        case EXPRINFO_POW:
          { double hh = __rval * __data.coeffs[eval_data.n];
            if(eval_data.n++ == 0)
              eval_data.r = hh+__data.params.nd();
            else
            { 
              if(hh == 0)
              {
                (*eval_data.d_data)[__data.node_num][0] = 0;
                (*eval_data.d_data)[__data.node_num][1] =
                                  log(eval_data.r)*__data.coeffs[1];
                eval_data.r = 1;
              }
              else
              {
                double h = pow(eval_data.r, hh);

                (*eval_data.d_data)[__data.node_num][0] =
                         hh*pow(eval_data.r, hh-1)*__data.coeffs[0];
                (*eval_data.d_data)[__data.node_num][1] =
                                       log(eval_data.r)*h*__data.coeffs[1];
                eval_data.r = h;
              }
            }
          }
          break;
        case EXPRINFO_EXP:
          { double h = exp(__rval*__data.coeffs[0]+__data.params.nd());
            eval_data.r = h;
            (*eval_data.d_data)[__data.node_num][0] = h*__data.coeffs[0];
          }
          break;
        case EXPRINFO_LOG:
          { double h = __rval*__data.coeffs[0]+__data.params.nd();
            eval_data.r = log(h);
            (*eval_data.d_data)[__data.node_num][0] = __data.coeffs[0]/h;
          }
          break;
        case EXPRINFO_SIN:
          { double h = __rval*__data.coeffs[0]+__data.params.nd();
            eval_data.r = sin(h);
            (*eval_data.d_data)[__data.node_num][0] = __data.coeffs[0]*cos(h);
          }
          break;
        case EXPRINFO_COS:
          { double h = __rval*__data.coeffs[0]+__data.params.nd();
            eval_data.r = cos(h);
            (*eval_data.d_data)[__data.node_num][0] = -__data.coeffs[0]*sin(h);
          }
          break;
        case EXPRINFO_ATAN2:
          { double hh = __rval * __data.coeffs[eval_data.n];
            if(eval_data.n++ == 0)
              eval_data.r = hh;
            else
            { double h = eval_data.r;
              h *= h;
              h += hh*hh;
              (*eval_data.d_data)[__data.node_num][0] = __data.coeffs[0]*hh/h;
              (*eval_data.d_data)[__data.node_num][1] =
                                                  -__data.coeffs[1]*eval_data.r/h;
              eval_data.r = atan2(eval_data.r,hh);
            }
          }
          break;
        case EXPRINFO_GAUSS:
          { double h = (__data.coeffs[0]*__rval-__data.params.d()[0])/
                                __data.params.d()[1];
            double k = exp(-h*h);
            eval_data.r = k;
            (*eval_data.d_data)[__data.node_num][0] =
                                 -2*__data.coeffs[0]*k*h/__data.params.d()[1];
          }
          break;
        case EXPRINFO_POLY:
          std::cerr << "func_d_evaluator: Polynomes NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_LIN:
        case EXPRINFO_QUAD:
          // we will never be 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;
            }
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          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
              (*eval_data.d_data)[__data.node_num][1] = 0.;
            }
            else
              (*eval_data.d_data)[__data.node_num][2] = 0.;
            (*eval_data.d_data)[__data.node_num][0] = 0.;
          }
          else // this is either the if or the else part
          {
            eval_data.r = __x;
            (*eval_data.d_data)[__data.node_num][eval_data.n] =
                                                __data.coeffs[eval_data.n];
            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
            }
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          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
            }
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          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;
            // for a piecewise constant function, 0 is always a subgradient
            (*eval_data.d_data)[__data.node_num][0] = 0.;
          }
          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
              }
              // for a piecewise constant function, 0 is always a subgradient
              (*eval_data.d_data)[__data.node_num][0] = 0.;
              (*eval_data.d_data)[__data.node_num][1] = 0.;
            }
            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;
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          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.u.p)->begin();
                _b != ((std::vector<double>*)eval_data.u.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.u.p)->push_back(__x);
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          eval_data.n++;
          if(eval_data.n == __data.n_children || ret == -1)
            delete (std::vector<double>*) eval_data.u.p;
          break;
        case EXPRINFO_HISTOGRAM:
          std::cerr << "func_d_evaluator: histogram NYI" << std::endl;
          throw "NYI";
          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;
            }
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          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;
              }
            }
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          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;
            }
          }
          // for a piecewise constant function, 0 is always a subgradient
          if(eval_data.n == 0)
          {
            std::vector<double>& v((*eval_data.d_data)[__data.node_num]);
            v.erase(v.begin(),v.end());
            v.insert(v.begin(),__data.n_children,0.);
          }
          eval_data.n++;
          break;
        case EXPRINFO_EXPECTATION:
          std::cerr << "func_d_evaluator: E NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_INTEGRAL:
          std::cerr << "func_d_evaluator: INT NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_LOOKUP:
        case EXPRINFO_PWLIN:
        case EXPRINFO_SPLINE:
        case EXPRINFO_PWCONSTLC:
        case EXPRINFO_PWCONSTRC:
          std::cerr << "func_d_evaluator: Table Operations NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_DET:
        case EXPRINFO_COND:
        case EXPRINFO_PSD:
        case EXPRINFO_MPROD:
        case EXPRINFO_FEM:
          std::cerr << "func_d_evaluator: Matrix Operations NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_RE:
        case EXPRINFO_IM:
        case EXPRINFO_ARG:
        case EXPRINFO_CPLXCONJ:
          std::cerr << "func_d_evaluator: Complex Numbers NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_CMPROD:
        case EXPRINFO_CGFEM:
          std::cerr << "func_d_evaluator: Const Matrix Operations NYI" << std::endl;
          throw "NYI";
          break;
        default:
          std::cerr << "func_d_evaluator: unknown function type " <<
                  __data.operator_type << std::endl;
          throw "Programming error";
          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,
                            &(*eval_data.d_data)[__data.node_num]);
    // use the cache only if it is worthwhile
    if(eval_data.f_cache && __data.n_parents > 1 && __data.n_children > 0)
      (*eval_data.f_cache)[__data.node_num] = eval_data.r;
    return ret;
  }

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

struct der_eval_type
{
  std::vector<std::vector<double> >* d_data;
  std::vector<std::vector<double> >* d_cache;
  std::vector<double>* grad_vec;
  const model* mod;
  double mult;
  double mult_trans;
  unsigned int child_n;
};

class der_eval : public
  cached_backward_evaluator_base<der_eval_type,expression_node,bool,
                                 model::walker>
{
private:
  typedef cached_backward_evaluator_base<der_eval_type,expression_node,
                                         bool,model::walker> _Base;

protected:
  bool is_cached(const node_data_type& __data)
     { 
       if(eval_data.d_cache && __data.n_parents > 1 && __data.n_children > 0
          && (*eval_data.d_cache)[__data.node_num].size() > 0 &&
          v_ind->match(__data.var_indicator()))
       {
         return true;
       }
       else
         return false;
     }

public:
  der_eval(std::vector<std::vector<double> >& __der_data, variable_indicator& __v,
           const model& __m, std::vector<std::vector<double > >* __d, 
           std::vector<double>& __grad)
  { // __grad MUST BE INITIALIZED WITH 0 AND THE CORRECT LENGTH
    eval_data.d_data = &__der_data;
    eval_data.d_cache = __d;
    eval_data.mod = &__m;
    eval_data.grad_vec = &__grad;
    eval_data.mult_trans = 1;
    eval_data.mult = 0;
    v_ind = &__v;
  }

  der_eval(const der_eval& __d) { eval_data = __d.eval_data; }

  ~der_eval() {}

  void new_point(std::vector<std::vector<double> >& __der_data,
                 const variable_indicator& __v)
  {
    eval_data.d_data = &__der_data;
    v_ind = &__v;
  }

  void new_result(std::vector<double>& __grad)
  {
    eval_data.grad_vec = &__grad;
  }

  void set_mult(double scal)
  {
    eval_data.mult_trans = scal;
  }
public:

  model::walker short_cut_to(const expression_node& __data)
      { return eval_data.mod->node(0); }  // return anything - not needed

  // at virtual nodes (ground, sky)
  void initialize()
  {
    eval_data.child_n = 0;
  }

  // before the children are visited
  int calculate(const expression_node& __data)
  {
    if(__data.operator_type == EXPRINFO_CONSTANT)
      return 0;
    else if(__data.operator_type == EXPRINFO_VARIABLE)
    {
      // update the gradient
      (*eval_data.grad_vec)[__data.params.nn()] += eval_data.mult_trans;
      return 0;
    }
    else if(__data.operator_type == EXPRINFO_LIN)
    { // update the gradient
      linalg_add(linalg_scale(eval_data.mod->lin[__data.params.nn()], eval_data.mult_trans),
          *eval_data.grad_vec,*eval_data.grad_vec);
      return 0;     // don't perform the remaining graph walk
    }
    else if(__data.operator_type == EXPRINFO_QUAD)
    { // update the gradient
      linalg_ssum(*eval_data.grad_vec, 2*eval_data.mult_trans,
                  (*eval_data.d_data)[__data.node_num]);
      return 0;     // don't perform the remaining graph walk
    }
    else if(__data.ev && (*__data.ev)[DER_EVALUATOR])
    // is there a short-cut evaluator defined?
    {
      linalg_ssum(*eval_data.grad_vec, eval_data.mult, 
          (*(der_evaluator)(*__data.ev)[DER_EVALUATOR])(
                              (*eval_data.d_data)[__data.node_num], *v_ind));
      return 0;     // don't perform the remaining graph walk
    }
    else if(eval_data.mult_trans == 0)
    // does not contribute to grad
      return 0;
    else
    {
      eval_data.child_n = 1;
      eval_data.mult = eval_data.mult_trans;
      if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.d_cache)
      { 
        eval_data.mult_trans = (*eval_data.d_data)[__data.node_num][0];
      }
      else // set multiplier for first child
        eval_data.mult_trans *= (*eval_data.d_data)[__data.node_num][0];
      return 1;
    }
  }

  // after all children have finished
  void cleanup(const expression_node& __data)
  {
    // use the cache only if it is worthwhile
    if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.d_cache
       && (*eval_data.d_cache)[__data.node_num].size() == 0)
    {
      (*eval_data.d_cache)[__data.node_num] = *eval_data.grad_vec;
      linalg_smult(*eval_data.grad_vec, eval_data.mult);
    }
  }

  void retrieve_from_cache(const expression_node& __data)
  {
    // the parallel f_cache MUST be initialized, too
    linalg_ssum(*eval_data.grad_vec, eval_data.mult_trans, 
          (*eval_data.d_cache)[__data.node_num]);
  }

  int update(const bool& __rval)
  {
    eval_data.child_n++;
    return 0;
  }

  // after the current child is processed
  int update(const expression_node& __data, const bool& __rval)
  {
    if(__data.n_children == 0)
      return 0;
    if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.d_cache)
    { 
      if(eval_data.child_n < __data.n_children)
        eval_data.mult_trans =
                (*eval_data.d_data)[__data.node_num][eval_data.child_n];
    }
    else if(eval_data.child_n < __data.n_children)
    { // prepare multiplier for the next child
      eval_data.mult_trans = eval_data.mult *
                (*eval_data.d_data)[__data.node_num][eval_data.child_n];
    }
    eval_data.child_n++;
    return 0;
  }

  bool calculate_value(bool eval_all)
  {
    return true;
  }
};

#endif /* _DER_EVALUATOR_H_ */

---