ihess_evaluator.h

Go to the documentation of this file.

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

// $Id: ihess_evaluator.h 679 2008-07-07 16:27:21Z 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 _IHESS_EVALUATOR_H_
#define _IHESS_EVALUATOR_H_

#define DEBUG_IHESS_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 <diffI.h>

namespace coco {

using namespace vgtl;

//                
// 1.) Preparation of the data-structures by the ihessPreparationEvaluator.
//      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.) ihessForwardEvaluator: 
//      fills the vector-pair (d,h) with the correct 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
//
// 3.) ihessBackwardEvaluator:
//      Here the stored values (D_{k,l}) and (M_{k,l}) are added up and multiplied
//      according to the recursion scheme.
//

struct ihessPreparationEvaluatorType
{
  std::vector<std::vector<interval> >* d; 
  std::vector<std::vector<interval> >* h; 
};

class  ihessPreparationEvaluator: public
  cached_forward_evaluator_base<ihessPreparationEvaluatorType, expression_node,
  bool, expression_const_walker>
{
private:
  typedef cached_forward_evaluator_base<ihessPreparationEvaluatorType,
                              expression_node,bool,expression_const_walker> _Base;

public:
  ihessPreparationEvaluator(std::vector<std::vector<interval> >& __d,
                            std::vector<std::vector<interval> >& __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<interval>());
     
  /*  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<interval>());
 */
    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<interval>());
 }

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

  ~ihessPreparationEvaluator() {}

  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, interval(0.));
    (*eval_data.h)[__data.node_num].insert(
        (*eval_data.h)[__data.node_num].end(), __data.n_children, interval(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 interval argument + preparation of derivative data
// for later backwards evaluation

struct ihessForwardEvaluatorReturnValue
{
  diffI f;
 
  unsigned int nn;
  tristate in_chn;
};

struct ihessForwardEvaluatorType
{
  const std::vector<diffI>* x;
  const std::vector<interval>* range;   
  std::vector<diffI>* f_cache;
  std::vector<std::vector<interval> >* h_data;
  std::vector<std::vector<interval> >* d_data; 
  std::vector<unsigned int>* t;
  std::vector<bool>* b;
  ihessForwardEvaluatorReturnValue r;
  diffI carry;
  const model* mod;
  void* p;
  interval d;
  unsigned int n, info; // needed for derivative of max, min
};

class ihessForwardEvaluator : public
  cached_forward_evaluator_base<ihessForwardEvaluatorType, expression_node,
                        ihessForwardEvaluatorReturnValue, expression_const_walker>
{
private:
  typedef cached_forward_evaluator_base<ihessForwardEvaluatorType,expression_node,
                 ihessForwardEvaluatorReturnValue, 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:
  diffI __power(double __coeff, diffI __x, int __exp)
  {
    if(__exp == 0)
      return 1.;
    else
    {
      diffI k = __coeff *__x;
      switch(__exp)
      {
        case 1:
          return k;
          break;
        case 2:
          return square(k);
          break;
        case -1:
          return 1./k;
          break;
        case -2:
          return 1./square(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(double 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:
  ihessForwardEvaluator(const std::vector<diffI>& __x,
              const std::vector<interval>& __rg,
              const variable_indicator& __v,
              const model& __m, std::vector<std::vector<interval> >& __d,
              std::vector<std::vector<interval> >& __h,
              std::vector<diffI>* __c) : _Base()
  {  
    #if DEBUG_IHESS_EVALUATOR
        cout<<"ihessForwardEvaluator.constructor"<< std::endl;
    #endif      
    eval_data.x = &__x;
    eval_data.range = &__rg;
    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 ihessPreparationEvaluator
    
    eval_data.h_data = &__h; // Note that eval_data.h_data is just the vector prepared by the ihessPreparationEvaluator
    //
    // what the forward-evaluator essentially does is to store the gradient-data resp. ihessian-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 ihessForwardEvaluatorReturnValue& __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.d = 0;
    v_ind = &__v;
     
  }

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

  ~ihessForwardEvaluator() {}

  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<diffI>& __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.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 newRange(const std::vector<interval>& __rg, const variable_indicator& __v)
  {
    eval_data.range = &__rg;
    v_ind = &__v;
  }

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

  int initialize(const expression_node& __data)
  {
    #if DEBUG_IHESS_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:
        case EXPRINFO_INVERT:
        case EXPRINFO_DIV:
        case EXPRINFO_CONSTANT:
          eval_data.r.f =diffI(__data.params.nd(),1);
          eval_data.carry=diffI(__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 = diffI(__data.params.nd(), 1);
          
          break;*/
        case EXPRINFO_IN:
        case EXPRINFO_AND:
        case EXPRINFO_NOGOOD:
          eval_data.r.f = diffI(1.,1);  
          break;
        case EXPRINFO_ALLDIFF:
       //   eval_data.u.p = (void*) new std::vector<interval>;
       //   ((std::vector<interval>*)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 = diffI(0.,1);
          break;
        case EXPRINFO_NORM:
          eval_data.info = 0;
          eval_data.r.f = diffI(0.,1);
          break;
        case EXPRINFO_DET:
        case EXPRINFO_PSD:
        // here construct square interval matrix
          break;
        case EXPRINFO_COND:
        case EXPRINFO_FEM:
        case EXPRINFO_MPROD:
        // here construct rectangular interval matrix
          break;
      }
      return 1;
    //}
  }

  void calculate(const expression_node& __data)
  {
         if(__data.operator_type > 0){
                 #if DEBUG_IHESS_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_IHESS_EVALUATOR
          std::cout<<"forward.retrieve_from_cache  node_num="<<eval_data.r.nn<< " node_counter="<< eval_data.n << std::endl;
          #endif
  }

  int update(const diffI& __rval)
  {
          #if DEBUG_IHESS_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 ihessForwardEvaluatorReturnValue& __rval) 
  {
          #if DEBUG_IHESS_EVALUATOR
          std::cout<<"forward.update  node_num="<<eval_data.r.nn<< " ";
          //std::cout<<"ihess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
          #endif
    int ret = 0;
    diffI __x(0.,1);
    diffI h(0.,1);
    diffI 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 = diffI(__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
          eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
          // 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.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            (*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_IHESS_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 diffI of the k-th input. 
            //
            //
            // h_data is calculated as in the der_evaluator.h
          
          if(eval_data.n == 0)
          {
            diffI 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
          {
            interval transfer_variable=eval_data.r.f[0];
            // interval 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;
            
            diffI 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++;
          if(eval_data.n == __data.n_children)
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
          break;
        }
       case EXPRINFO_MONOME:
       { 
               std::cerr << "ihessForwardEvaluator: Polynomes NYI" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_QUAD";
          break;
        }
        case EXPRINFO_MAX:
        {
          std::cerr << "ihessForwardEvaluator: MAX" << std::endl;
          throw api_exception(apiee_other,"Regularity error: ihessForwardEvaluator. EXPRINFO_MAX");
          break;
        }
        case EXPRINFO_MIN:
        { 
          std::cerr << "ihessForwardEvaluator: MIN" << std::endl;
          throw api_exception(apiee_other,"Regularity error: ihessForwardEvaluator. EXPRINFO_MIN");
          break;
        }
         
        case EXPRINFO_SCPROD:
        { 
           std::cerr << "ihessForwardEvaluator: SCPROD" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_SCPROD";
            break;
        }
        case EXPRINFO_NORM:
        { 
               std::cerr << "ihessForwardEvaluator: NORM" << std::endl;
            throw "NYI: ihessForwardEvaluator. 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_IHESS_EVALUATOR
                std::cout<<"EXPRINFO_INVERT"<<std::endl;
           #endif  
            
            h= 1 / __rval.f ;
            eval_data.r.f = eval_data.r.f * h; 
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            
            (*eval_data.d_data)[__data.node_num][0] = 
                       -__data.params.nd()*sqr(h[0]);
                       
            (*eval_data.h_data)[__data.node_num][0] = 
                       2*__data.params.nd()*sqr(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_IHESS_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();
            diffI h=  1/__rval.f;
            interval g1=  eval_data.r.f[1]; 
            interval temp= sqr(h[0]);
            eval_data.r.f = constant * eval_data.r.f *  h ;
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
          //   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 ihessian_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 diffI of the first input and (y[0],y[1]) is the diffI 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] = temp*(- constant * g1 + 2 * eval_data.r.f[0] *__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 = square(h);
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            #if DEBUG_IHESS_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_IHESS_EVALUATOR
                std::cout<<"EXPRINFO_INTPOWER"<<std::endl;
          #endif
            int hl = __data.params.nn();
            if(hl == 0)
            {
              eval_data.r.f = 1;
              eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
              
              (*eval_data.h_data)[__data.node_num][0] = 0;
              (*eval_data.d_data)[__data.node_num][0] = 0;
                                
            }
            else
            {
              diffI kl = __data.coeffs[0]*__rval.f;
              switch(hl)
              {
                case 1:
                  eval_data.r.f = kl;
                  eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
                  (*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 = square(kl);
                  eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
                        
                  (*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:
                  { diffI h = 1/kl;
                    h[0].intersectwith((*eval_data.range)[__data.node_num]);
                    eval_data.r.f = h;
                   (*eval_data.h_data)[__data.node_num][0] = 2*ipow(h[0],3)*kl[1]*__data.coeffs[0];
                   (*eval_data.d_data)[__data.node_num][0] = -sqr(h[0])*__data.coeffs[0];
                  }
                  break;
                case -2:
                  { 
                    diffI h = 1/kl;
                    diffI k = square(h);
                    k[0].intersectwith((*eval_data.range)[__data.node_num]);
                    eval_data.r.f = k;
                 
                    (*eval_data.h_data)[__data.node_num][0] =  6*sqr(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:
                  {

                    diffI k=__power(1.,kl, hl-1);
                    eval_data.r.f = k*kl;
                    eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
                    
                    (*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_IHESS_EVALUATOR
                std::cout<<"EXPRINFO_SQROOT"<<std::endl;
          #endif
            h = sqrt(__data.coeffs[0]*__rval.f+__data.params.nd());
            h[0].intersectwith((*eval_data.range)[__data.node_num]);
            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:
        { 
          std::cerr << "ihessForwardEvaluator: ABS" << std::endl;
          throw api_exception(apiee_other,"Regularity error: ihessForwardEvaluator. EXPRINFO_ABS");
          break;
        }
        case EXPRINFO_POW: 
        { 
          //std::cerr << "ihessForwardEvaluator: EXPRINFO_POW noch nicht implementiert, da 'diffI pow(const diffI &a, const diffI &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
          diffI 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)
              {
                      
                interval 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
              {*/
            interval x =  eval_data.r.f[0];
            interval logx = log(eval_data.r.f[0]);
            interval y = __rval.f[0] * __data.coeffs[1];
            interval yy1;  // y(y-1)
            if(y.contains(0.5)) {
              interval u = abs(y-0.5);
              u = u.sup();
              u *= u-1;
              yy1 = interval(-0.25, u.sup());
            } else {
              interval u = abs(y-0.5);
              yy1 = u.inf();
              yy1 *= yy1-1;
              u = u.sup();
              u *= u-1;
              yy1.hullwith(u);
            }
            //diffI h = std::pow(eval_data.r.f, hh);
            eval_data.r.f = power(eval_data.r.f, __rval.f * __data.coeffs[1]);
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            interval temp2 = eval_data.r.f[0] / sqr(x);
                
            (*eval_data.h_data)[__data.node_num][0] =
               temp2 * ( yy1 * __data.coeffs[0] * __data.coeffs[0] * eval_data.carry[1]
                         + x * (1+ y * logx ) * __data.coeffs[0] * __data.coeffs[1] * __rval.f[1] ) ;
            temp2 = eval_data.r.f[0] / x;
            (*eval_data.h_data)[__data.node_num][1] =
              temp2 *  ( (1+ y * logx ) * __data.coeffs[0] * __data.coeffs[1] * eval_data.carry[1]
                         + sqr(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
            diffI h = exp(__rval.f*__data.coeffs[0]+__data.params.nd());
            h[0].intersectwith((*eval_data.range)[__data.node_num]);
            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
            diffI h = __rval.f*__data.coeffs[0]+__data.params.nd();
            interval k=__data.coeffs[0]/h[0];
            eval_data.r.f = log(h);
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            (*eval_data.d_data)[__data.node_num][0]  = k;
            (*eval_data.h_data)[__data.node_num][0] = - sqr(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        
            diffI h = __rval.f*__data.coeffs[0]+__data.params.nd();
            eval_data.r.f = sin(h);
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            interval cosh0 = sqrt(1-sqr(eval_data.r.f[0]));
            cosh0.intersectwith(cos(h[0]));
           
            (*eval_data.d_data)[__data.node_num][0] =__data.coeffs[0] * cosh0;
            (*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          
            diffI h = __rval.f*__data.coeffs[0]+__data.params.nd();
            eval_data.r.f = cos(h);
            eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            interval sinh0 = sqrt(1-sqr(eval_data.r.f[0]));
            sinh0.intersectwith(sin(h[0]));
        
            (*eval_data.d_data)[__data.node_num][0] = -__data.coeffs[0]*sinh0;
            (*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          
            diffI hh = __rval.f * __data.coeffs[eval_data.n];
            if(eval_data.n++ == 0) {
              eval_data.carry=__rval.f;   
              eval_data.r.f = hh;
            }
            else
            { 
                interval h = eval_data.r.f[0];
                h = sqr(h);
                h = h +  sqr(hh[0]);
                interval gradx = hh[0]/h;
                interval 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; 
                interval ihessdiag= 2 * gradx * grady;
                // interval ihessoffdiag= __data.coeffs[0]*__data.coeffs[1] *(gradx+grady)*(gradx-grady);
                interval ihessoffdiag= __data.coeffs[0]*(__data.coeffs[1]*(sqr(hh[0])-sqr(eval_data.r.f[0]))/sqr(h));
                (*eval_data.h_data)[__data.node_num][0] =  __data.coeffs[0]* __data.coeffs[0] * ihessdiag * eval_data.carry[1] 
                                                                +   ihessoffdiag * __rval.f[1];
                (*eval_data.h_data)[__data.node_num][1] =   ihessoffdiag * eval_data.carry[1]
                                                                 -__data.coeffs[1]* __data.coeffs[1] * ihessdiag * __rval.f[1] ;
                eval_data.r.f = atan2(eval_data.r.f,hh);
                eval_data.r.f[0].intersectwith((*eval_data.range)[__data.node_num]);
            }
          }
          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        
                  
                  diffI h = (__data.coeffs[0]*__rval.f-__data.params.d()[0])/
                                __data.params.d()[1];
            diffI k = exp(-h*h);
            k[0].intersectwith((*eval_data.range)[__data.node_num]);
            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*sqr(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 << "ihessForwardEvaluator: POLY" << std::endl;
            throw "NYI: ihessForwardEvaluator. 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 << "ihessForwardEvaluator: QUAD" << std::endl;
          throw "NYI: ihessForwardEvaluator. EXPRINFO_QUAD";
          break;
        }
        case EXPRINFO_IN:
          {
            std::cerr << "ihessForwardEvaluator: IN" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_IN";
            break;
          }
        case EXPRINFO_IF:
         {
            std::cerr << "ihessForwardEvaluator: IF" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_IF";
            break;
          }
        case EXPRINFO_AND:
         {
            std::cerr << "ihessForwardEvaluator: AND" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_AND";
            break;
          }
        case EXPRINFO_OR:
         {
            std::cerr << "ihessForwardEvaluator: OR" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_OR";
            break;
          }
        case EXPRINFO_NOT:
           {
            std::cerr << "ihessForwardEvaluator: NOT" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_NOT";
            break;
          }
          
        case EXPRINFO_IMPLIES:
         {
            std::cerr << "ihessForwardEvaluator: IMPLIES" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_IMPLIES";
            break;
          }
        case EXPRINFO_COUNT:
          {
            std::cerr << "ihessForwardEvaluator: COUNT" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_COUNT";
            break;
          }
        case EXPRINFO_ALLDIFF:
          {
            std::cerr << "ihessForwardEvaluator: ALLDIFF" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_ALLDIFF";
            break;
          }
        case EXPRINFO_HISTOGRAM:
         {
            std::cerr << "ihessForwardEvaluator: HISTOGRAM" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_HISTOGRAM";
            break;
          }
        case EXPRINFO_LEVEL:
          {
            std::cerr << "ihessForwardEvaluator: LEVEL" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_LEVEL";
            break;
          }
        case EXPRINFO_NEIGHBOR:
          {
            std::cerr << "ihessForwardEvaluator: NEIGHBOR" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_NEIGHBOR";
            break;
          }
        case EXPRINFO_NOGOOD:
          {
            std::cerr << "ihessForwardEvaluator: NOGOOD" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_NOGOOD";
            break;
          }
        case EXPRINFO_EXPECTATION:
          {
            std::cerr << "ihessForwardEvaluator: EXPECTATION" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_EXPECTATION";
            break;
          }
        case EXPRINFO_INTEGRAL:
          {
            std::cerr << "ihessForwardEvaluator: INTEGRAL" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_INTEGRAL";
            break;
          }
        case EXPRINFO_LOOKUP:
        case EXPRINFO_PWLIN:
        case EXPRINFO_SPLINE:
        case EXPRINFO_PWCONSTLC:
        case EXPRINFO_PWCONSTRC:
          {
            std::cerr << "ihessForwardEvaluator:  SPLINE_STUFF" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_SPLINE_STUFF";
            break;
          }
        case EXPRINFO_DET:
        case EXPRINFO_COND:
        case EXPRINFO_PSD:
        case EXPRINFO_MPROD:
        case EXPRINFO_FEM:
          {
            std::cerr << "ihessForwardEvaluator: MATRIX_STUFF" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_MATRIX_STUFF";
            break;
          }
        case EXPRINFO_RE:
        case EXPRINFO_IM:
        case EXPRINFO_ARG:
        case EXPRINFO_CPLXCONJ:
          {
            std::cerr << "ihessForwardEvaluator: COMPLEX_STUFF" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_COMPLEX_STUFF";
            break;
          }
        case EXPRINFO_CMPROD:
        case EXPRINFO_CGFEM:
          {
            std::cerr << "ihessForwardEvaluator: FEM_STUFF" << std::endl;
            throw "NYI: ihessForwardEvaluator. EXPRINFO_FEM_STUFF";
            break;
          }
        default:
        {
          std::cerr << "ihessForwardEvaluator: 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;
  }

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

 


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

class  ihessBackwardEvaluator : public
  cached_backward_evaluator_base<ihessBackwardEvaluatorType,expression_node,bool,
                                 expression_const_walker>
{
private:
  typedef cached_backward_evaluator_base<ihessBackwardEvaluatorType,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:
  ihessBackwardEvaluator(std::vector<std::vector<interval> >*__der_data, 
            std::vector<std::vector<interval> >* __hess_data,
            variable_indicator& __v,
            const model& __m, 
            std::vector<std::vector<interval > >* __d, 
            std::vector<std::vector<interval > >* __h, 
            std::vector<interval>* __grad,
            std::vector<interval>* __hess)
  { // __grad MUST BE INITIALIZED WITH 0 AND THE CORRECT LENGTH
    #if DEBUG_IHESS_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;
  }
  
 
  ihessBackwardEvaluator(const ihessBackwardEvaluator& __he) { eval_data = __he.eval_data; }

  ~ihessBackwardEvaluator() {}

  void newPoint(std::vector<std::vector<interval> >* __der_data,
                 std::vector<std::vector<interval> >* __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<interval>* __grad, std::vector<interval>* __hess)
  {
     eval_data.grad_vec = __grad;
     eval_data.hess_vec = __hess;
  }

  void set_mult(interval 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_IHESS_EVALUATOR
          std::cout<<"backward.initialize" << std::endl;
          //std::cout<<"ihess_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_IHESS_EVALUATOR
                        std::cout<<"backward.calculate:  EXPRINFO_CONSTANT" << std::endl;
                                //std::cout<<"ihess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
                  #endif  
                  break;
                  
                  }
                  case EXPRINFO_VARIABLE:
                  {
                  #if DEBUG_IHESS_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_IHESS_EVALUATOR
                        std::cout<<"backward.calculate:  EXPRINFO_LIN" << std::endl;
                                //std::cout<<"ihess_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_IHESS_EVALUATOR
                        std::cout<<"backward.calculate:  EXPRINFO_QUAD" << std::endl;
                                //std::cout<<"ihess_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_IHESS_EVALUATOR
                        
                                std::cout<<"backward.calculate:  DEFAULT optype="<<__data.operator_type << std::endl;
                #endif  
            
                eval_data.child_n = 1;
                

                        #if DEBUG_IHESS_EVALUATOR
                                //std::cout<<"oldvalue hmult=" << eval_data.h_mult << std::endl;
                                //std::cout<<"oldvalue hmultrans=" <<  eval_data.h_mult_trans << std::endl;
                                //std::cout<<"ihess_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_IHESS_EVALUATOR
                                //std::cout<<"newvalue hmult=" << eval_data.h_mult << std::endl;
                                //std::cout<<"newvalue hmultrans=" <<  eval_data.h_mult_trans << std::endl;
                                //std::cout<<"ihess_evaluator.update vecsize="<< (*eval_data.h_data)[__data.node_num].size() << std::endl;
                        #endif  
                        
                
                
                
                if(eval_data.d_mult_trans != 0. ) {
                        #if DEBUG_IHESS_EVALUATOR
                                //std::cout<<"oldvalue dmult=" << eval_data.d_mult << std::endl;
                                //std::cout<<"oldvalue dmultrans=" <<  eval_data.d_mult_trans << std::endl;
                                //std::cout<<"ihess_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_IHESS_EVALUATOR
                                //std::cout<<"newvalue dmult=" << eval_data.d_mult << std::endl;
                                //std::cout<<"newvalue dmultrans=" <<  eval_data.d_mult_trans << std::endl;
                                //std::cout<<"ihess_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_IHESS_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_IHESS_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_IHESS_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_IHESS_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_IHESS_EVALUATOR
                std::cout<<"backward.update"<<std::endl;;
          #endif  
          if(__data.n_children == 0) {
                  #if DEBUG_IHESS_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_IHESS_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_IHESS_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_IHESS_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_IHESS_EVALUATOR
                std::cout<<"backward.calculate_value" << std::endl;
          #endif  
    return true;
  }
};



} // namespace coco

#endif /* _IHESS_EVALUATOR_H_ */

---