hess_evaluator.h

Go to the documentation of this file.

// Hessian Evaluator implementation -*- C++ -*-

// $Id: hess_evaluator.h 732 2008-10-17 11:33:07Z mihaly $
// Copyright (C) 2007-2008 Simon Haller, 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 _HESS_EVALUATOR_H_
#define _HESS_EVALUATOR_H_

#define DEBUG_HESS_EVALUATOR 0
#define DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 0

#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>
#include <diffNumber.h>

namespace coco {

using namespace vgtl;

//                
// here a little introduction, to hessematrix-vector multiplication evaluator:
// a short sketch of the algorithm can be found in "../docs/extra/hessian_vector_product_evaluator.dvi"
// note that all the variable-names are taken from "../api/der_evaluator.h", so in order to understand
// the "hess_evaluator.h" source, it is maybe better to start with simpler "der_eavluator.h". 
// Nevertheless at least i've added some comments..
// From the structural point of view both evaluators are almost identical. 
//
// the hessian vector multiplication is carried out in 3 steps:
//
// 1.) Preparation of the data-structures by the hessPreparationEvaluator. Once for each function-tree only!
//      For each node a pair of vectors (d,h) is created of the size of its inputs
//      Note that d[j] corresponds to the (D_{k,l}) in the dvi-file describing the algorithm (each j can be described by a pair k,l), while
//      h[j] corresponds to (M_{k,l}), so both vectors are of size n_{k-1} (by notation used in the dvi-file) which is the number of "inputs" for node j.        
//              
// 2.) hessForwardEvaluator: 
//      fills the vector-pair (d,h) with the right data for all nodes. Here the complicated part of the algotithm happens. 
//      But actually we just have to calculate (D_{k,l}) and (M_{k,l}) for each node (k,l) <-> j
//      according to the formulas given in the dvi-file!
//
// 3.) hessBackwardEvaluator:
//      Here the stored values (D_{k,l}) and (M_{k,l}) are added up and multiplied according to recursion scheme derived in
//      in the last part of the dvi-file! Check out the formulas for the "P_{k,l} and Q_{k,l}"
//

struct hessPreparationEvaluatorType
{
  std::vector<std::vector<double> >* d; 
  std::vector<std::vector<double> >* h; 
};

class  hessPreparationEvaluator: public
  cached_forward_evaluator_base<hessPreparationEvaluatorType, expression_node, bool, expression_const_walker>
{
private:
  typedef cached_forward_evaluator_base<hessPreparationEvaluatorType,
                                    expression_node,bool,expression_const_walker> _Base;

public:
  hessPreparationEvaluator(std::vector<std::vector<double> >& __d, std::vector<std::vector<double> >& __h,
              unsigned int _num_of_nodes)
  {
    eval_data.d=&__d;
    eval_data.h=&__h;
    if((*eval_data.d).size() < _num_of_nodes && eval_data.d!=NULL)
      (*eval_data.d).insert((*eval_data.d).end(), _num_of_nodes-(*eval_data.d).size(), std::vector<double>());
     
  /*  if((*eval_data.h).size() < _num_of_nodes)
      (*eval_data.h).insert((*eval_data.h).end(),
                      _num_of_nodes-(*eval_data.h).size(), linalg::matrix<double>());
 */
    if((*eval_data.h).size() < _num_of_nodes  )
       (*eval_data.h).insert((*eval_data.h).end(), _num_of_nodes-(*eval_data.h).size(), std::vector<double>());
 }

  hessPreparationEvaluator(const  hessPreparationEvaluator& __x) { eval_data.d = __x.eval_data.d;
                                        eval_data.h = __x.eval_data.h; }

  ~hessPreparationEvaluator() {}

  void initialize() { return; }

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

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

  int initialize(const expression_node& __data)
  {
    
    (*eval_data.d)[__data.node_num].insert(
        (*eval_data.d)[__data.node_num].end(), __data.n_children, 0);
    (*eval_data.h)[__data.node_num].insert(
        (*eval_data.h)[__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; }
        
 /* void clearAll() 
  {
          
  }*/

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

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

struct hessForwardEvaluatorReturnValue
{
  diffNumber f;
 
  unsigned int nn;
  tristate in_chn;
};

struct hessForwardEvaluatorType
{
  const std::vector<diffNumber>* x;
  std::vector<diffNumber>* f_cache;
  std::vector<std::vector<double> >* h_data;
  std::vector<std::vector<double> >* d_data; 
  std::vector<unsigned int>* t;
  std::vector<bool>* b;
  hessForwardEvaluatorReturnValue r;
  diffNumber carry;
  const model* mod;
  union { void* p; double d; } u;
  unsigned int n, 
           info; // needed for derivative of max, min
};

class hessForwardEvaluator : public
  cached_forward_evaluator_base<hessForwardEvaluatorType, expression_node,
                                hessForwardEvaluatorReturnValue, expression_const_walker>
{
private:
  typedef cached_forward_evaluator_base<hessForwardEvaluatorType,expression_node,
                                        hessForwardEvaluatorReturnValue, expression_const_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:
  diffNumber __power(double __coeff, diffNumber __x, int __exp)
  {
    if(__exp == 0)
      return 1.;
    else
    {
      diffNumber 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:
          return power(k,__exp);
         // 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;
      }
    }
    return 0; // just to make the compiler happy
  }

  void __calc_max(diffNumber h, const expression_node& __data)
  {
    if(h[0] >= eval_data.r.f[0])
    {
      (*eval_data.d_data)[__data.node_num][eval_data.info] = 0;
      if(h[0] > eval_data.r.f[0])
        (*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.f = h;
      eval_data.info = eval_data.n;
    }
    else
    {
      (*eval_data.d_data)[__data.node_num][eval_data.n] = 0;
    }
    (*eval_data.h_data)[__data.node_num][eval_data.n] = 0;
  }

public:
  hessForwardEvaluator(const std::vector<diffNumber>& __x, const variable_indicator& __v,
              const model& __m, std::vector<std::vector<double> >& __d,
              std::vector<std::vector<double> >& __h,
              std::vector<diffNumber>* __c) : _Base()
  {  
    #if DEBUG_HESS_EVALUATOR
        cout<<"hessForwardEvaluator.constructor"<< std::endl;
    #endif      
    eval_data.x = &__x;
    eval_data.f_cache = __c;
    eval_data.mod = &__m;
    eval_data.d_data = &__d; 
    // Note that eval_data.d_data is just the vector prepared by the hessPreparationEvaluator
    
    eval_data.h_data = &__h; // Note that eval_data.h_data is just the vector prepared by the hessPreparationEvaluator
    //
    // what the forward-evaluator essentially does is to store the gradient-data resp. hessian-data at a specific node (nr. j) of the graph into the
    // d__data[j] resp. h_data[j];
    // Note that d__data[j] corresponds to the (D_{k,l}) in the dvi-file describing the algorithm (each j can be described by a pair k,l), while
    // h__data[j] corresponds to (M_{k,l}), so both vectors are of size n_{k-1} (by notation used in the dvi-file) which is the number of "inputs" for node j.
    //
    //
    // just go to the function 
    // "int update(const expression_node& __data, const hessForwardEvaluatorReturnValue& __rval)"
    // to see where the funny things happens :). Essentially 
    // "case EXPRINFO_DIV:" is good a place to start with, as it has two "inputs" which means that it is not too trivial!
    // for further comments proceed there.... 
    //
    
    eval_data.n = 0;
    eval_data.r.f.setDeg(1);
    eval_data.carry.setDeg(1);
    eval_data.u.d = 0;
    v_ind = &__v;
     
  }

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

  ~hessForwardEvaluator() {}

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

  void newPoint(const std::vector<diffNumber>& __x, const variable_indicator& __v)
  {
    eval_data.x = &__x;
  //  std::cout << "new point: " << __x << std::endl;
    v_ind = &__v;
    eval_data.n = 0;
    eval_data.r.f[0]=0.;eval_data.r.f[1]=0.; 
    eval_data.carry[0]=0.; eval_data.carry[1]=0.;
    eval_data.u.d = 0;
    v_ind = &__v;

    
#if 0
    for(unsigned int i=0; i<(*eval_data.d_data).size(); i++)
    {
            for(unsigned int j=0; j< ((*eval_data.d_data)[i]).size(); j++)
                    (*eval_data.d_data)[i][j] =0.0; 
    }
    for(unsigned int i=0; i<(*eval_data.h_data).size(); i++)
            for(unsigned int j=0; j< ((*eval_data.h_data)[i]).size(); j++)
                    (*eval_data.h_data)[i][j] =0.0;
#endif
    
  }

  void initialize() {
          #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"forward.initialize" << std::endl;
          #endif
          #if DEBUG_HESS_EVALUATOR
                        std::cout<<"forward.initialize" << std::endl;
                #endif  
  return; }

  int initialize(const expression_node& __data)
  {
    #if DEBUG_HESS_EVALUATOR
        std::cout<<"forward.initialize node_num="<<eval_data.r.nn<< " node_counter="<< eval_data.n << std::endl;
    #endif      
    eval_data.r.nn = __data.node_num;
    eval_data.r.in_chn = __data.sem.info_flags.has_0chnbase;
    eval_data.n = 0;
 
    switch(__data.operator_type)
      {
        case EXPRINFO_VARIABLE:
        
           eval_data.r.f = (*eval_data.x)[__data.params.nn()];
          break;
        case EXPRINFO_SUM:
        case EXPRINFO_PROD:
        case EXPRINFO_MAX:
        case EXPRINFO_MIN:
          eval_data.info = 0;
          // no break == continue
        case EXPRINFO_INVERT:
        case EXPRINFO_DIV:
        case EXPRINFO_CONSTANT:
          eval_data.r.f =diffNumber(__data.params.nd(),1);
          eval_data.carry=diffNumber(__data.params.nd(),1);
          break;      
              
        /*case EXPRINFO_MAX:
        case EXPRINFO_MIN:
          eval_data.info = 0;
          // no break == continue
        case EXPRINFO_SUM:
        case EXPRINFO_PROD:
        case EXPRINFO_INVERT:
        case EXPRINFO_DIV:      
          std::cout<< "init="<<__data.params.nd()<<";" ;
          eval_data.r.f = diffNumber(__data.params.nd(), 1);
          
          break;*/
        case EXPRINFO_IN:
        case EXPRINFO_AND:
        case EXPRINFO_NOGOOD:
          eval_data.r.f = diffNumber(1.,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.f = diffNumber(0.,1);
          break;
        case EXPRINFO_NORM:
          eval_data.info = 0;
          eval_data.r.f = diffNumber(0.,1);
          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){
                 #if DEBUG_HESS_EVALUATOR
                        std::cout<<"forward.calculate node_num="<<eval_data.r.nn<< " node_counter="<< eval_data.n << std::endl;
                #endif
         }  
          
  
  }

  void retrieve_from_cache(const expression_node& __data)
  {
          #if DEBUG_HESS_EVALUATOR
          std::cout<<"forward.retrieve_from_cache  node_num="<<eval_data.r.nn<< " node_counter="<< eval_data.n << std::endl;
          #endif
  }

  int update(const diffNumber& __rval)
  {
          #if DEBUG_HESS_EVALUATOR
          std::cout<<"forward.update  node_num="<<eval_data.r.nn<< " node_counter="<< eval_data.n << std::endl;
          #endif
    eval_data.r.f = __rval;
    return 0;
  }

  
  int update(const expression_node& __data, const hessForwardEvaluatorReturnValue& __rval) 
  {
          #if DEBUG_HESS_EVALUATOR
          std::cout<<"forward.update  node_num="<<eval_data.r.nn<< " ";
          //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
          #endif
    int ret = 0;
    diffNumber __x(0.,1);
    diffNumber h(0.,1);
    diffNumber k(0.,1);
    //if(eval_data.r.in_chn == t_false && __rval.in_chn != t_false)
    //  (*eval_data.t).push_back(__rval.nn);
    
        
    if(__data.operator_type < 0)
    {
      switch(__data.operator_type)
      {
        case EXPRINFO_CONSTANT:
        {
                
          //(*eval_data.d_data)[__data.node_num][eval_data.n]  = 0.;
          //(*eval_data.d_data)[__data.node_num][eval_data.n]  = 0.;
          eval_data.r.f =diffNumber( __data.params.nd(),1);
          #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_CONSTANT (__data.params.nd()="<< __data.params.nd() << "), ";
          #endif
           
         
          // don't care about d_data, CONSTANTS are ALWAYS leafs
        break;
      }
        case EXPRINFO_VARIABLE:
        {
         
        //  (*eval_data.d_data)[__data.node_num][eval_data.n]  = 1.;
        //  (*eval_data.d_data)[__data.node_num][eval_data.n]  = 0.;
        #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_VARIABLE (__data.params.nn()="<< __data.params.nn() << "), ";
          #endif
         
          // don't care about d_data, VARIABLES are ALWAYS leafs
        break;
      }
        case EXPRINFO_SUM:
        
        case EXPRINFO_MEAN:{
        {
                double hd = __data.coeffs[eval_data.n];
            #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_MEAN (__data.params.coeffs="<< __data.coeffs[eval_data.n] << ";__data.params.hd()="<< __data.params.nd() <<"), ";
          #endif
        
            eval_data.r.f = eval_data.r.f + hd * __rval.f;             // this is just the standard function evaluation
            (*eval_data.d_data)[__data.node_num][eval_data.n]  = hd;   // here we keep track of the [__data.node_num]
           
            //for(int i = 0 ; i < eval_data.n; i++)
                    
            (*eval_data.h_data)[__data.node_num][eval_data.n] = 0.; 
            
             #if DEBUG_HESS_EVALUATOR
                std::cout<<"EXPRINFO_MEAN " <<  eval_data.r.f << std::endl;
                //std::cout<<"EXPRINFO_MEAN:" << eval_data.r.f << ","<< hd<< ","<< __rval.f <<std::endl;
                #endif 
            eval_data.n++;
        }
        break;
      }
        case EXPRINFO_PROD: {
        
        
          #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_PROD (__data.params.nd()="<< __data.params.nd() << ";__data.coeffs[0]="<< __data.coeffs[0] <<"), ";
          #endif
          
          //
            // This one is a bit complicated! But what essentially happens is:
            // h_data is calculated by using the formula HESSIAN ( of the node (x1,....,xn) -> x1 * x2 * ... * xn ) x DIRECTIONAL DERIVATIVE 
            // ( 0                          , x3[0] * ... * xn[0]              , x2[0]* x4[0] * ... xn[0]               , .....   , x2[0] * x3[0] * ... * x(n-1)[0]         )  *  ( x1[1] )
            // ( x3[0] * ... * xn[0]        , 0                                , x1[0]* x4[0] * ... xn[0]               , .....   , x1[0] * x3[0] * ... * x(n-1)[0]         )  *  ( x2[1] )     
            // ( x2[0] * x4[0] *... * xn[0] , x1[0] * x4[0] * ... * xn[0]      , 0                                      , .....   , x1[0] * x2[0] * x4[0] * ... * x(n-1)[0] )  *  ( x3[1] )
            // (      ................................................................................................................................                      )  *  (.......)
            // ( x2[0] * ... * x(n-1)[0]    ,  x1[0] * x3[0] * ... * x(n-1)[0] ,  x1[0] * x2[0] * x4[0] ... * x(n-1)[0] , ......  , 0                                       )  *  ( xn[1] )
            // where (xk[0],xk[1],) is the diffNumber of the k-th input. 
            //
            //
            // h_data is calculated as in the der_evaluator.h
          
          if(eval_data.n == 0)
          {
            diffNumber temp(1);
            temp=eval_data.r.f * __rval.f;
            eval_data.r.f = temp;
            
            eval_data.carry =__rval.f;
           // (*eval_data.h_data)[__data.node_num][0] = __data.params.nd();
            //if(eval_data.d_data!=NULL)
            //std::cout << "paramprod: " << temp <<std::endl;
            (*eval_data.d_data)[__data.node_num][0] = __data.params.nd();
            (*eval_data.h_data)[__data.node_num][0] = 0;
          }
          else
          {
            double transfer_variable=eval_data.r.f[0];
           // double old_ddata= 
          (*eval_data.d_data)[__data.node_num][eval_data.n] = transfer_variable;
          (*eval_data.h_data)[__data.node_num][eval_data.n] = 
                (*eval_data.h_data)[__data.node_num][eval_data.n-1]* eval_data.carry[0]  + (*eval_data.d_data)[__data.node_num][eval_data.n-1] * eval_data.carry[1]; 
           
           
            //(*eval_data.h_data)[__data.node_num][eval_data.n] =0;
            //for(int i=0; i<eval_data.n; i++)
            //(*eval_data.h_data)[__data.node_num][eval_data.n] = eval_data.r.f[1];
            //std::cout << "&&& "  << eval_data.carry << " " << __rval.f << std::endl;
            
            eval_data.carry= __rval.f;
            
            diffNumber temp(1);
            temp=eval_data.r.f * __rval.f;
            eval_data.r.f = temp;
            //for(int i = eval_data.n-1; i >= 0; i--)
            //    (*eval_data.h_data)[__data.node_num][i] 
            //          = eval_data.r.f[1]* __rval.f[1];
             //if(eval_data.d_data!=NULL)  
             
             for(int i = eval_data.n-1; i >= 0; i--) {
                (*eval_data.h_data)[__data.node_num][i] 
                                = (*eval_data.h_data)[__data.node_num][i] * __rval.f[0];
                                
                 
                (*eval_data.h_data)[__data.node_num][i] 
                        = (*eval_data.h_data)[__data.node_num][i] + (*eval_data.d_data)[__data.node_num][i]  *  __rval.f[1]  ;  
                        
                (*eval_data.d_data)[__data.node_num][i] 
                        = (*eval_data.d_data)[__data.node_num][i]  * __rval.f[0]; // * __rval.f[0];
                //for(int j=0; j < i; i++) transfer_variable *=
                
                
                }
                
                
                
             }
              
         
          eval_data.n++;
          break;
        }
       case EXPRINFO_MONOME:
       { 
               std::cerr << "hessForwardEvaluator: Polynomes NYI" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_QUAD";
          break;
        }
        case EXPRINFO_MAX:
        {
          __calc_max(__rval.f * __data.coeffs[eval_data.n], __data);
          ++eval_data.n;
          break;
        }
        case EXPRINFO_MIN:
        { 
          { diffNumber h = __rval.f * __data.coeffs[eval_data.n];
            if(h[0] <= eval_data.r.f[0])
            {
              (*eval_data.d_data)[__data.node_num][eval_data.info] = 0;
              if(h[0] < eval_data.r.f[0])
                (*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.f = h;
              eval_data.info = eval_data.n;
            }
            else
            {
              (*eval_data.d_data)[__data.node_num][eval_data.n] = 0;
            }
            (*eval_data.h_data)[__data.node_num][eval_data.n] = 0;
          }
          ++eval_data.n;
          break;
        }
         
        case EXPRINFO_SCPROD:
        { 
           std::cerr << "hessForwardEvaluator: SCPROD" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_SCPROD";
            break;
        }
        case EXPRINFO_NORM:
        { 
               std::cerr << "hessForwardEvaluator: NORM" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_NORM";
          break;
        }
        case EXPRINFO_INVERT:
        {
         #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_INVERT (__data.params.nd()="<< __data.params.nd() << "), ";
          #endif
          
          #if DEBUG_HESS_EVALUATOR
                std::cout<<"EXPRINFO_INVERT"<<std::endl;
           #endif  
            
            h= 1 / __rval.f ;
            eval_data.r.f =eval_data.r.f * h; 
            
                    (*eval_data.d_data)[__data.node_num][0] = 
                       -__data.params.nd()*h[0]*h[0];
                       
                    (*eval_data.h_data)[__data.node_num][0] = 
                       2*__data.params.nd()*h[0]* h[0] * h[0] * __rval.f[1];
          
        break;
        }
        case EXPRINFO_DIV:  
        { // has two "inputs"
                #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_DIV (__data.params.nd()="<< __data.params.nd() <<";__data.coeffs[0]="<<__data.coeffs[0]<< "), ";
          #endif
        #if DEBUG_HESS_EVALUATOR
                std::cout<<"EXPRINFO_DIV"<<std::endl;
          #endif  
          if(eval_data.n++ == 0) {
                  eval_data.r.f =  __rval.f;
          }
          else
          {
            double constant=__data.params.nd();
            diffNumber h=  1/__rval.f;
            double g1=  eval_data.r.f[1]; 
            double temp=h[0] * h[0];
            eval_data.r.f = constant * eval_data.r.f *  h ;
          //   std::cout << constant << std::endl;
             /* 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;   der-evaluator   */
            //
            // compare these expressions with the "M_{k,l}" from the hessian_vector_product_evaluator.dvi file
            // It calculated by the formula HESSIAN ( of the node (x,y) -> x/y ) x DIRECTIONAL DERIVATIVE 
            // (  0                 -1/y[0]² ) * ( x[1] )
            // (-1/y[0]²         2x[0]/y[0]² )   ( y[1] )
            // where (x[0],x[1]) is the diffNumber of the first input and (y[0],y[1]) is the diffNumber of the second input. 
            //

            (*eval_data.h_data)[__data.node_num][0] = - constant * temp *  __rval.f[1]  ;
            (*eval_data.h_data)[__data.node_num][1] = - constant * g1 * temp + 2 * eval_data.r.f[0] * temp *__rval.f[1];
            
            (*eval_data.d_data)[__data.node_num][0] = h[0] * constant ;
            (*eval_data.d_data)[__data.node_num][1] = - eval_data.r.f[0] * h[0];
            
          }
  
        break;
        }
        case EXPRINFO_SQUARE:
        {
            #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_SQUARE (__data.params.nd()="<< __data.params.nd() <<";__data.coeffs[0]=" << __data.coeffs[0] << "), ";
            #endif
            h = __data.coeffs[0]*__rval.f+__data.params.nd();
            eval_data.r.f = h*h;
             #if DEBUG_HESS_EVALUATOR
                std::cout<<"EXPRINFO_SQUARE " << eval_data.r.f  <<std::endl;
             #endif
           
            (*eval_data.h_data)[__data.node_num][0] = 2 * h[1] * __data.coeffs[0];
            (*eval_data.d_data)[__data.node_num][0] =  2 * h[0] * __data.coeffs[0];
            
          break;
        }
        case EXPRINFO_INTPOWER:
        {
           #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_INTPOWER (__data.params.nn()="<< __data.params.nn() << "; __data.coeffs="<<__data.coeffs<<"), ";
           #endif
           #if DEBUG_HESS_EVALUATOR
                std::cout<<"EXPRINFO_INTPOWER"<<std::endl;
          #endif
            int hl = __data.params.nn();
            if(hl == 0)
            {
              eval_data.r.f = 1;
              
              (*eval_data.h_data)[__data.node_num][0] = 0;
              (*eval_data.d_data)[__data.node_num][0] = 0;
                                
            }
            else
            {
              diffNumber kl = __data.coeffs[0]*__rval.f;
              switch(hl)
              {
                case 1:
                  eval_data.r.f = kl;
                   (*eval_data.h_data)[__data.node_num][0] = 0;
                   (*eval_data.d_data)[__data.node_num][0] = __data.coeffs[0];
                  break;
                case 2:
                  eval_data.r.f = kl*kl;
                 
                        
                  (*eval_data.h_data)[__data.node_num][0] =  2*kl[1]*__data.coeffs[0];
                 
                  (*eval_data.d_data)[__data.node_num][0] = 2*kl[0]*__data.coeffs[0];
                        
                  break;
                case -1:
                  { diffNumber h = 1/kl;
                    eval_data.r.f = h;
                   (*eval_data.h_data)[__data.node_num][0] = 2*h[0]*h[0]*h[0]*kl[1]*__data.coeffs[0];
                   (*eval_data.d_data)[__data.node_num][0] = -h[0]*h[0]*__data.coeffs[0];       
                  }
                  break;
                case -2:
                  { 
                    diffNumber h = 1/kl;
                    diffNumber k = h*h;
                    eval_data.r.f = k;
                 
                    (*eval_data.h_data)[__data.node_num][0] =  6*k[0]*k[0]*kl[1]*__data.coeffs[0];
                  
                     (*eval_data.d_data)[__data.node_num][0] = -2*h[0]*k[0]*__data.coeffs[0];
                   
                  }
                  break;
                default:
                  {

                    diffNumber k=__power(1.,kl, hl-1);
                    eval_data.r.f = k*kl;
                    
                    (*eval_data.h_data)[__data.node_num][0] = hl * __data.coeffs[0] * k[1];
                    (*eval_data.d_data)[__data.node_num][0] = hl * __data.coeffs[0] * k[0];
                  }
                  break;
              }
            }
          
          break;
        }
        case EXPRINFO_SQROOT:
        {
            #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_SQROOT(__data.params.nd()="<< __data.params.nd() << "; __data.coeffs[0]="<<__data.coeffs<<"), ";
           #endif
           #if DEBUG_HESS_EVALUATOR
                std::cout<<"EXPRINFO_SQROOT"<<std::endl;
          #endif
            h = sqrt(__data.coeffs[0]*__rval.f+__data.params.nd());
            eval_data.r.f = h;
            k = 0.5 * __data.coeffs[0] / h;
            
             (*eval_data.h_data)[__data.node_num][0] =  k[1];
                    if(eval_data.d_data!=NULL) {
                           (*eval_data.d_data)[__data.node_num][0] =  k[0];
                    }
          
          break;
        }
        case EXPRINFO_ABS:
        { 
          diffNumber hd = __data.coeffs[0] * __rval.f;
        
          if(hd[0] < 0)
            eval_data.r.f = -hd;
          else
            eval_data.r.f = hd;
          (*eval_data.d_data)[__data.node_num][0] =
                    hd[0] > 0 ? __data.coeffs[0] : (hd[0] < 0 ? -__data.coeffs[0] : 0);
          
          (*eval_data.h_data)[__data.node_num][0] = 0.; 
            
          #if DEBUG_HESS_EVALUATOR
            std::cout<<"EXPRINFO_ABS " <<  eval_data.r.f << std::endl;
            //std::cout<<"EXPRINFO_MEAN:" << eval_data.r.f << ","<< hd<< ","<< __rval.f <<std::endl;
          #endif 
        }
        case EXPRINFO_POW: 
        { 
          //std::cerr << "HessForwardEvaluator: EXPRINFO_POW noch nicht implementiert, da 'diffNumber pow(const diffNumber &a, const diffNumber &b)' nicht existiert." << std::endl;
          //throw "NYI";
          //break;
                 
          #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
            std::cout<<"EXPRINFO_POW(n="<<eval_data.n<<";__data.params.nd()="<< __data.params.nd() << "; __data.coeffs[0]="<<__data.coeffs<<"), ";
          #endif
          diffNumber hh = __rval.f * __data.coeffs[eval_data.n];
          if(eval_data.n++ == 0)
          {
            eval_data.carry = __rval.f;
            eval_data.r.f = hh+__data.params.nd();
          }
          else
          { 
              /*if(hh == 0)
              {
                      
                double temp= std::log(eval_data.r.f)*__data.coeffs[1];
                
                (*eval_data.d_data)[__data.node_num][0] = 0;
                (*eval_data.d_data)[__data.node_num][1] = temp;
                               
                (*eval_data.h_data)[__data.node_num][0]= __data.coeffs[1] * __data.coeffs[1] *  
                (*eval_data.h_data)[__data.node_num][1]= __datacoeffs[0] * __datacoeffs[1]/eval_data.r.f *__rval.carry[1] + temp * temp * __rval.f[1];            
                eval_data.r.f = 1;
              }
              else
              {*/
            double x =   eval_data.r.f[0];
            double logx =    std::log(eval_data.r.f[0]);
            double y =    __rval.f[0] * __data.coeffs[1];
            //diffNumber h = std::pow(eval_data.r.f, hh);
            eval_data.r.f = power(eval_data.r.f, __rval.f * __data.coeffs[1]);
            double temp2=eval_data.r.f[0] / (x*x);
                
            (*eval_data.h_data)[__data.node_num][0] =
               temp2 * ( y * (y-1) * __data.coeffs[0] * __data.coeffs[0] * eval_data.carry[1]
                         + x * (1+ y * logx ) * __data.coeffs[0] * __data.coeffs[1] * __rval.f[1] ) ;
            temp2 *= x;
            (*eval_data.h_data)[__data.node_num][1] =
              temp2 *  ( (1+ y * logx ) * __data.coeffs[0] * __data.coeffs[1] * eval_data.carry[1]
                         + logx * logx * x * __data.coeffs[1] * __data.coeffs[1] * __rval.f[1]);
                
            (*eval_data.d_data)[__data.node_num][0] = y *  temp2; 
                         
            (*eval_data.d_data)[__data.node_num][1] = logx * eval_data.r.f[0] * __data.coeffs[1];
                         
          }
          break;
        }
        case EXPRINFO_EXP:
          {
            #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                std::cout<<"EXPRINFO_EXP(__data.params.nd()="<< __data.params.nd() << "; __data.coeffs[0]="<<__data.coeffs[0]<<"), ";
            #endif
            diffNumber h = exp(__rval.f*__data.coeffs[0]+__data.params.nd());
            eval_data.r.f = h;
            (*eval_data.d_data)[__data.node_num][0] = __data.coeffs[0]* h[0];
            (*eval_data.h_data)[__data.node_num][0] = __data.coeffs[0] *__data.coeffs[0] * h[0] * __rval.f[1];
                
            
          
          break;
          }
        case EXPRINFO_LOG:
          {  
                #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                        std::cout<<"EXPRINFO_LOG(__data.params.nd()="<< __data.params.nd() << "; __data.coeffs[0]="<<__data.coeffs[0]<<"), ";
                #endif
                  diffNumber h = __rval.f*__data.coeffs[0]+__data.params.nd();
                  double k=__data.coeffs[0]/h[0];
                  eval_data.r.f = log(h);
                  (*eval_data.d_data)[__data.node_num][0]  = k;
                   (*eval_data.h_data)[__data.node_num][0] = - k * k * __rval.f[1];
          }
          break;
       case EXPRINFO_SIN:
          { 
            #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                        std::cout<<"EXPRINFO_SIN(__data.params.nd()="<< __data.params.nd() << "; __data.coeffs[0]="<<__data.coeffs[0]<<"), ";
            #endif        
            diffNumber h = __rval.f*__data.coeffs[0]+__data.params.nd();
            eval_data.r.f = sin(h);
           
            (*eval_data.d_data)[__data.node_num][0] =__data.coeffs[0] * std::cos(h[0]);
            (*eval_data.h_data)[__data.node_num][0] = - __data.coeffs[0] * __data.coeffs[0] * eval_data.r.f[0] * __rval.f[1];
          }
          break;
        case EXPRINFO_COS:
          { 
            #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                        std::cout<<"EXPRINFO_COS(__data.params.nd()="<< __data.params.nd() << "; __data.coeffs[0]="<<__data.coeffs[0]<<"), ";
            #endif          
            diffNumber h = __rval.f*__data.coeffs[0]+__data.params.nd();
            eval_data.r.f = cos(h);
        
            (*eval_data.d_data)[__data.node_num][0] = -  __data.coeffs[0] * std::sin(h[0]);
            (*eval_data.h_data)[__data.node_num][0] = - __data.coeffs[0]* __data.coeffs[0] * eval_data.r.f[0] *__rval.f[1];
          }
          break;
         case EXPRINFO_ATAN2:
          { 
                #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                        std::cout<<"EXPRINFO_ATAN2(__data.params.nd()="<< __data.params.nd() << "; __data.coeffs[0]="<<__data.coeffs[0]<<"), ";
                #endif      
                
                diffNumber hh = __rval.f * __data.coeffs[eval_data.n];
                if(eval_data.n++ == 0) {
              eval_data.carry=__rval.f;   
                eval_data.r.f = hh;}
            else
            { 
                double h = eval_data.r.f[0];
                h = h * h;
                h = h +  hh[0]*hh[0];
                double gradx = hh[0]/h;
                double grady =  eval_data.r.f[0]/h;
                
                (*eval_data.d_data)[__data.node_num][0] = __data.coeffs[0]   * gradx;
                (*eval_data.d_data)[__data.node_num][1] = - __data.coeffs[1] * grady; 
                double hessdiag= 2 * gradx * grady;
                double hessoffdiag= __data.coeffs[0]* __data.coeffs[1] *(gradx+grady)*(gradx-grady);
                (*eval_data.h_data)[__data.node_num][0] =  __data.coeffs[0]* __data.coeffs[0] * hessdiag * eval_data.carry[1] 
                                                                +   hessoffdiag * __rval.f[1];
                (*eval_data.h_data)[__data.node_num][1] =   hessoffdiag * eval_data.carry[1]
                                                                 -__data.coeffs[1]* __data.coeffs[1] * hessdiag * __rval.f[1] ;
                eval_data.r.f = atan2(eval_data.r.f,hh);
            }
          }
          break;
        case EXPRINFO_GAUSS:
          { 
                #if DEBUG_FHE_LIST_NODES_WITH_PARAMETERS 
                        std::cout<<"EXPRINFO_GAUSS(__data.params.d()="<< __data.params.d() << "; __data.coeffs[0]="<<__data.coeffs[0]<<"), ";
                #endif        
                  
                  diffNumber h = (__data.coeffs[0]*__rval.f-__data.params.d()[0])/
                                __data.params.d()[1];
            diffNumber k = exp(-h*h);
            eval_data.r.f = k;
           
            (*eval_data.d_data)[__data.node_num][0] =
                                 -2*h[0]*__data.coeffs[0]*k[0]/__data.params.d()[1];
                (*eval_data.h_data)[__data.node_num][0] =
                                 2*(2*h[0]*h[0]-2) * __data.coeffs[0] *__data.coeffs[0]  * k[0] / (__data.params.d()[1]*__data.params.d()[1]);           
          }
          break;
          
        case EXPRINFO_POLY:
         {
            std::cerr << "hessForwardEvaluator: POLY" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_POLY";
            break;
          }
        /*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 be here!
          std::cerr << "hessForwardEvaluator: QUAD" << std::endl;
          throw "NYI: hessForwardEvaluator. EXPRINFO_QUAD";
          break;
        }
        case EXPRINFO_IN:
          {
            std::cerr << "hessForwardEvaluator: IN" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_IN";
            break;
          }
        case EXPRINFO_IF:
         {
            std::cerr << "hessForwardEvaluator: IF" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_IF";
            break;
          }
        case EXPRINFO_AND:
         {
            std::cerr << "hessForwardEvaluator: AND" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_AND";
            break;
          }
        case EXPRINFO_OR:
         {
            std::cerr << "hessForwardEvaluator: OR" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_OR";
            break;
          }
        case EXPRINFO_NOT:
           {
            std::cerr << "hessForwardEvaluator: NOT" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_NOT";
            break;
          }
          
        case EXPRINFO_IMPLIES:
         {
            std::cerr << "hessForwardEvaluator: IMPLIES" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_IMPLIES";
            break;
          }
        case EXPRINFO_COUNT:
          {
            std::cerr << "hessForwardEvaluator: COUNT" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_COUNT";
            break;
          }
        case EXPRINFO_ALLDIFF:
          {
            std::cerr << "hessForwardEvaluator: ALLDIFF" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_ALLDIFF";
            break;
          }
        case EXPRINFO_HISTOGRAM:
         {
            std::cerr << "hessForwardEvaluator: HISTOGRAM" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_HISTOGRAM";
            break;
          }
        case EXPRINFO_LEVEL:
          {
            std::cerr << "hessForwardEvaluator: LEVEL" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_LEVEL";
            break;
          }
        case EXPRINFO_NEIGHBOR:
          {
            std::cerr << "hessForwardEvaluator: NEIGHBOR" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_NEIGHBOR";
            break;
          }
        case EXPRINFO_NOGOOD:
          {
            std::cerr << "hessForwardEvaluator: NOGOOD" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_NOGOOD";
            break;
          }
        case EXPRINFO_EXPECTATION:
          {
            std::cerr << "hessForwardEvaluator: EXPECTATION" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_EXPECTATION";
            break;
          }
        case EXPRINFO_INTEGRAL:
          {
            std::cerr << "hessForwardEvaluator: INTEGRAL" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_INTEGRAL";
            break;
          }
        case EXPRINFO_LOOKUP:
        case EXPRINFO_PWLIN:
        case EXPRINFO_SPLINE:
        case EXPRINFO_PWCONSTLC:
        case EXPRINFO_PWCONSTRC:
          {
            std::cerr << "hessForwardEvaluator:  SPLINE_STUFF" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_SPLINE_STUFF";
            break;
          }
        case EXPRINFO_DET:
        case EXPRINFO_COND:
        case EXPRINFO_PSD:
        case EXPRINFO_MPROD:
        case EXPRINFO_FEM:
          {
            std::cerr << "hessForwardEvaluator: MATRIX_STUFF" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_MATRIX_STUFF";
            break;
          }
        case EXPRINFO_RE:
        case EXPRINFO_IM:
        case EXPRINFO_ARG:
        case EXPRINFO_CPLXCONJ:
          {
            std::cerr << "hessForwardEvaluator: COMPLEX_STUFF" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_COMPLEX_STUFF";
            break;
          }
        case EXPRINFO_CMPROD:
        case EXPRINFO_CGFEM:
          {
            std::cerr << "hessForwardEvaluator: FEM_STUFF" << std::endl;
            throw "NYI: hessForwardEvaluator. EXPRINFO_FEM_STUFF";
            break;
          }
        default:
        {
          std::cerr << "HessForwardEvaluator: 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.f = __data.f_evaluate(eval_data.n++, __data.params.nn(),
                        *eval_data.x, *v_ind, eval_data.r.f, __rval.f,
                            &(*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.f;
    return ret;
  }

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

 


struct hessBackwardEvaluatorType
{
  std::vector<std::vector<double> >* d_data;
  std::vector<std::vector<double> >* h_data;
  std::vector<std::vector<double> >* d_cache;
  std::vector<std::vector<double> >* h_cache;
  std::vector<double>* grad_vec;
  std::vector<double>* hess_vec;
  bool calc_grad;
  const model* mod;
  double d_mult;
  double hessian_of_current_node; // out of order
  bool hess_zero;
  double d_mult_trans;
  double h_mult;
  double h_mult_trans;
  bool is_linear;
  unsigned int child_n;
  //unsigned int h_child_n;
  //unsigned int d_child_n;
};

class  hessBackwardEvaluator : public
  cached_backward_evaluator_base<hessBackwardEvaluatorType,expression_node,bool,
                                 expression_const_walker>
{
private:
  typedef cached_backward_evaluator_base<hessBackwardEvaluatorType,expression_node,
                                         bool,expression_const_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:
  hessBackwardEvaluator(std::vector<std::vector<double> >*__der_data, 
            std::vector<std::vector<double> >* __hess_data,
            variable_indicator& __v,
            const model& __m, 
            std::vector<std::vector<double > >* __d, 
            std::vector<std::vector<double > >* __h, 
            std::vector<double>* __grad,
            std::vector<double>* __hess)
  { // __grad MUST BE INITIALIZED WITH 0 AND THE CORRECT LENGTH
    #if DEBUG_HESS_EVALUATOR
                        cout<<"backward.constructor"<< std::endl;
    #endif      
    eval_data.d_data = __der_data;
    eval_data.h_data = __hess_data;
    eval_data.d_cache = __d;
    eval_data.h_cache = __h;
    eval_data.mod = &__m;
    eval_data.grad_vec = __grad;
    eval_data.hess_vec = __hess;
    
    eval_data.d_mult_trans = 1.;
    eval_data.d_mult = 0.;
    eval_data.h_mult_trans = 0.;
    eval_data.h_mult = 0.;
    eval_data.hess_zero =true;
    eval_data.hessian_of_current_node =0.;
    v_ind = &__v;
  }
  
 
  hessBackwardEvaluator(const hessBackwardEvaluator& __he) { eval_data = __he.eval_data; }

  ~hessBackwardEvaluator() {}

  void newPoint(std::vector<std::vector<double> >* __der_data,
                 std::vector<std::vector<double> >* __hess_data,
                 const variable_indicator& __v)
  {
    if(__der_data!=NULL) eval_data.d_data = __der_data;
    else eval_data.calc_grad=false;
    eval_data.h_data = __hess_data;
    v_ind = &__v;
  }

  void newResult(std::vector<double>* __grad, std::vector<double>* __hess)
  {
     eval_data.grad_vec = __grad;
     eval_data.hess_vec = __hess;
  }

  void set_mult(double scal)
  {
    eval_data.d_mult_trans = scal;
    eval_data.h_mult_trans = 0.;
  }
public:

  expression_const_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()
  {
        
          #if DEBUG_HESS_EVALUATOR
          std::cout<<"backward.initialize" << std::endl;
          //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
          #endif
    eval_data.child_n = 0;
  }

  // before the children are visited
  int calculate(const expression_node& __data)
  {
        
          switch(__data.operator_type) { 
                  case EXPRINFO_CONSTANT: 
                  {
                  #if DEBUG_HESS_EVALUATOR
                        std::cout<<"backward.calculate:  EXPRINFO_CONSTANT" << std::endl;
                                //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
                  #endif  
                  break;
                  
                  }
                  case EXPRINFO_VARIABLE:
                  {
                  #if DEBUG_HESS_EVALUATOR
                        std::cout<<"backward.calculate:  EXPRINFO_VARIABLE" << std::endl;
                  #endif  
     
                        int i=__data.params.nn();
                        
                        (*eval_data.grad_vec)[i] = (*eval_data.grad_vec)[i]+ eval_data.d_mult_trans;
                        (*eval_data.hess_vec)[i] = (*eval_data.hess_vec)[i] +eval_data.h_mult_trans;
                        eval_data.hessian_of_current_node=0.;
                        
                
                  break;
                  }
        case EXPRINFO_SUM: 
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_MEAN:
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_PROD:
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_MONOME:
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_MAX:
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_SCPROD:
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_NORM:
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_INVERT:
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_DIV:  // has two children
        {eval_data.hessian_of_current_node=0.;}
        case EXPRINFO_SQUARE:
        {eval_data.hessian_of_current_node=2.;}
    
    /*else if(__data.operator_type == EXPRINFO_LIN)
    { // update the gradient
             #if DEBUG_HESS_EVALUATOR
                        std::cout<<"backward.calculate:  EXPRINFO_LIN" << std::endl;
                                //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
            #endif  
              linalg::linalg_add(linalg_scale(eval_data.mod->lin[__data.params.nn()], eval_data.d_mult_trans),
                        *eval_data.grad_vec,*eval_data.grad_vec);
                 
                linalg::linalg_add(linalg_scale(eval_data.mod->lin[__data.params.nn()], eval_data.h_mult_trans),
                        *eval_data.hess_vec,*eval_data.hess_vec);
                        
      return 0;     // don't perform the remaining graph walk
    }**/
   /* else if(__data.operator_type == EXPRINFO_QUAD)
    { // update the gradient
                  #if DEBUG_HESS_EVALUATOR
                        std::cout<<"backward.calculate:  EXPRINFO_QUAD" << std::endl;
                                //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
            #endif  
                  linalg::linalg_ssum(*eval_data.grad_vec, 2*eval_data.d_mult_trans,
                                        (*eval_data.d_data)[__data.node_num]);
                  linalg::linalg_ssum(*eval_data.hess_vec, 2*eval_data.h_mult_trans,
                                        (*eval_data.h_data)[__data.node_num]);
      return 0;     // don't perform the remaining graph walk
    }*/
    default: 
    {
                #if DEBUG_HESS_EVALUATOR
                        
                                std::cout<<"backward.calculate:  DEFAULT optype="<<__data.operator_type << std::endl;
                #endif  
            
                eval_data.child_n = 1;
                

                        #if DEBUG_HESS_EVALUATOR
                                //std::cout<<"oldvalue hmult=" << eval_data.h_mult << std::endl;
                                //std::cout<<"oldvalue hmultrans=" <<  eval_data.h_mult_trans << std::endl;
                                //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
                        #endif  
                        
                        eval_data.h_mult = eval_data.h_mult_trans;
                        
                        eval_data.h_mult_trans  = eval_data.h_mult_trans * (*eval_data.d_data)[__data.node_num][0] 
                                                 +eval_data.d_mult_trans * (*eval_data.h_data)[__data.node_num][0]; 
                        // here we use the recursion derived in the last lines of  hessian_vector_product.dvi!
                        
                        
                        
                        //(*eval_data.h_data)[__data.node_num][0];
                        
                        #if DEBUG_HESS_EVALUATOR
                                //std::cout<<"newvalue hmult=" << eval_data.h_mult << std::endl;
                                //std::cout<<"newvalue hmultrans=" <<  eval_data.h_mult_trans << std::endl;
                                //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
                        #endif  
                        
                
                
                
                if(eval_data.d_mult_trans != 0. ) {
                        #if DEBUG_HESS_EVALUATOR
                                //std::cout<<"oldvalue dmult=" << eval_data.d_mult << std::endl;
                                //std::cout<<"oldvalue dmultrans=" <<  eval_data.d_mult_trans << std::endl;
                                //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
                        #endif  
                        
                        eval_data.d_mult = eval_data.d_mult_trans;
                        //eval_data.h_mult = eval_data.h_mult_trans;
                        eval_data.d_mult_trans *= (*eval_data.d_data)[__data.node_num][0];
    
                        #if DEBUG_HESS_EVALUATOR
                                //std::cout<<"newvalue dmult=" << eval_data.d_mult << std::endl;
                                //std::cout<<"newvalue dmultrans=" <<  eval_data.d_mult_trans << std::endl;
                                //std::cout<<"hess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
                        #endif  
        
                }
                 
    }
            
        
    }
          
        return 1;  
  }
  // after all children have finished
  void cleanup(const expression_node& __data) 
  {
           #if DEBUG_HESS_EVALUATOR
                std::cout<<"backward.cleanup" << std::endl;
          #endif
    // use the cache only if it is worthwhile
  /*   if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.h_cache
       && (*eval_data.h_cache)[__data.node_num].size() == 0)
    {
        #if DEBUG_HESS_EVALUATOR
                std::cout<<"backward.cleanup h weirdo cache op" << std::endl;
          #endif
     // (*eval_data.h_cache)[__data.node_num] = *eval_data.hess_vec;
    //  linalg::linalg_smult(*eval_data.hess_vec, eval_data.h_mult);
      // linalg::linalg_ssum( linalg::linalg_smult(*eval_data.hess_vec, eval_data.d_mult), *eval_data.grad_vec ,   eval_data.h_mult) ;
    }
    
    if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.d_cache
       && (*eval_data.d_cache)[__data.node_num].size() == 0)
    {
             #if DEBUG_HESS_EVALUATOR
                std::cout<<"backward.cleanup d weirdo cache op" << std::endl;
          #endif
         (*eval_data.d_cache)[__data.node_num] = *eval_data.grad_vec;
      
      linalg::linalg_smult(*eval_data.grad_vec, eval_data.d_mult);
    } */
  }

/*  void retrieve_from_cache(const expression_node& __data)
  {
    // the parallel f_cache MUST be initialized, too
    if(eval_data.grad_vec!=NULL) linalg::linalg_ssum(*eval_data.grad_vec, eval_data.d_mult_trans, 
          (*eval_data.d_cache)[__data.node_num]);
    linalg::linalg_ssum(*eval_data.hess_vec, eval_data.h_mult_trans, 
          (*eval_data.h_cache)[__data.node_num]);
  }*/

  int update(const bool& __rval)
  {
                   #if DEBUG_HESS_EVALUATOR
                std::cout<<"backward.update bool" << std::endl;
          #endif  
    eval_data.child_n++;
    return 0;
  }

  // after the current child is processed
  int update(const expression_node& __data, const bool& __rval) //hier muss ein fehler sein!
  {
             #if DEBUG_HESS_EVALUATOR
                std::cout<<"backward.update"<<std::endl;;
          #endif  
          if(__data.n_children == 0) {
                  #if DEBUG_HESS_EVALUATOR
                                //std::cout << "no children " <<  eval_data.h_mult_trans << std::endl;
                        #endif  
                         return 0;
          }
     
    
                /*if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.h_cache)
                { 
                        if(eval_data.child_n < __data.n_children) {
                        eval_data.h_mult_trans =  (*eval_data.h_data)[__data.node_num][eval_data.child_n];}
                         #if DEBUG_HESS_EVALUATOR
                                std::cout << "#1 " <<  eval_data.h_mult_trans << std::endl;
                        #endif  
                }
                else */
                if(eval_data.child_n < __data.n_children)
                { // prepare multiplier for the next child
                        
                         #if DEBUG_HESS_EVALUATOR
                                //std::cout<<""<< std::endl;
                                //std::cout << "\t ddata=" <<  (*eval_data.d_data)[__data.node_num][eval_data.child_n] << " hdata="<< (*eval_data.h_data)[__data.node_num][eval_data.child_n] <<std::endl;
                                //std::cout << "\t dmult=" << eval_data.d_mult  <<" hmult="<<   eval_data.h_mult  <<std::endl;
                                //std::cout << "\t dmulttrans=" << eval_data.d_mult_trans  <<" hmulttrans="<<   eval_data.h_mult_trans  <<std::endl;
                                //std::cout<<""<< std::endl;
                        #endif  
                        
                        eval_data.d_mult_trans = eval_data.d_mult * (*eval_data.d_data)[__data.node_num][eval_data.child_n];
                        
                        eval_data.h_mult_trans = eval_data.h_mult * (*eval_data.d_data)[__data.node_num][eval_data.child_n] 
                                                +eval_data.d_mult * (*eval_data.h_data)[__data.node_num][eval_data.child_n];
                        //eval_data.h_mult_trans=(*eval_data.h_data)[__data.node_num][eval_data.child_n];
                        #if DEBUG_HESS_EVALUATOR
                                //std::cout<<""<< std::endl;
                                //std::cout << "\t ddata=" <<  (*eval_data.d_data)[__data.node_num][eval_data.child_n] << " hdata="<< (*eval_data.h_data)[__data.node_num][eval_data.child_n] <<std::endl;
                                //std::cout << "\t dmult=" << eval_data.d_mult  <<" hmult="<<   eval_data.h_mult  <<std::endl;
                                //std::cout << "\t dmulttrans=" << eval_data.d_mult_trans  <<" hmulttrans="<<   eval_data.h_mult_trans  <<std::endl;
                                //std::cout<<""<< std::endl;
                        #endif  
                }               
                
        
    
    
    eval_data.child_n++;
    return 0;
  }

  bool calculate_value(bool eval_all)
  {
             #if DEBUG_HESS_EVALUATOR
                std::cout<<"backward.calculate_value" << std::endl;
          #endif  
    return true;
  }
};



} // namespace coco

#endif /* _HESS_EVALUATOR_H_ */

---