---

islp_evaluator.h

Go to the documentation of this file.

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

// Copyright (C) 2001-2003 Hermann Schichl
//
// This file is part of the COCONUT API.  This library
// is free software; you can redistribute it and/or modify it under the
// terms of the Library GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// Library GNU General Public License for more details.

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

#ifndef _ISLP_EVALUATOR_H_
#define _ISLP_EVALUATOR_H_

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

using namespace vgtl;

struct func_islp_return_type
{
  double f;
  interval rg;
  interval fi;
};

typedef bool (*prep_islp_evaluator)();
typedef func_islp_return_type (*func_islp_evaluator)(
                  const std::vector<interval>* __x, const variable_indicator& __v,
                  std::vector<interval>& __islp_data);
typedef std::vector<interval>& (*islp_evaluator)(
                const std::vector<interval>& __d_dat,
                const variable_indicator& __v);

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

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

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

  ~prep_islp_eval() {}

  void initialize() { return; }

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

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

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

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

  int update(bool __rval) { return 0; }

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

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

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

struct func_islp_eval_type
{
  const std::vector<double>* z;
  const std::vector<interval>* range;
  std::vector<double>* f;
  std::vector<std::vector<interval> >* islp_data;
  const model* mod;
  union { void *p; interval_st d; unsigned int info; } u;
  func_islp_return_type r;
  unsigned int n; 
};

class func_islp_eval : public
  cached_forward_evaluator_base<func_islp_eval_type, expression_node,
                                func_islp_return_type, model::walker>
{
private:
  typedef cached_forward_evaluator_base<func_islp_eval_type,expression_node,
                                   func_islp_return_type, model::walker> _Base;

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

public:
  func_islp_eval(const std::vector<double>& __z,
                 const std::vector<interval>& __rg,
                 const variable_indicator& __v, const model& __m,
                 std::vector<std::vector<interval> >& __d,
                 std::vector<double>& __f) : _Base()
  {
    eval_data.z = &__z;
    eval_data.range = &__rg;
    eval_data.f = &__f;
    eval_data.islp_data = &__d;
    eval_data.mod = &__m;
    eval_data.n = 0;
    eval_data.r.f = 0;
    eval_data.r.rg = 0;
    eval_data.r.fi = 0;
    v_ind = &__v;
  }

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

  ~func_islp_eval() {}

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

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

  void new_range(const std::vector<interval>& __rg,
                 const variable_indicator& __v)
  {
    eval_data.range = &__rg;
    v_ind = &__v;
  }

  void initialize() { return; }

  int initialize(const expression_node& __data)
  {
    eval_data.n = 0;
    if(__data.ev != NULL && (*__data.ev)[FUNC_ISLP_EVALUATOR] != NULL)
    // is there a short-cut evaluator defined?
    {
      eval_data.r = (*(func_islp_evaluator)(*__data.ev)[FUNC_ISLP_EVALUATOR])
                        (eval_data.range, *v_ind,
                         (*eval_data.islp_data)[__data.node_num]);
      return 0;     // don't perform the remaining graph walk
    }
    else
    {
      switch(__data.operator_type)
      {
        case EXPRINFO_MAX:
        case EXPRINFO_MIN:
          eval_data.u.info = 0;
          // no break == continue
        case EXPRINFO_SUM:
        case EXPRINFO_INVERT:
          eval_data.r.fi = eval_data.r.f = __data.params.nd();
          break;
        case EXPRINFO_PROD:
          eval_data.r.rg = eval_data.r.fi = eval_data.r.f = __data.params.nd();
          break;
        case EXPRINFO_IN:
        case EXPRINFO_AND:
        case EXPRINFO_NOGOOD:
          eval_data.r.fi = eval_data.r.f = 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:
          eval_data.r.fi = eval_data.r.f = 0.;
          break;
        case EXPRINFO_NORM:
          eval_data.r.fi = eval_data.r.f = 0.;
          eval_data.u.info = 0;
          break;
        case EXPRINFO_LEVEL:
          eval_data.u.p = (void*) new std::vector<interval>;
          ((std::vector<interval>*)eval_data.u.p)->reserve(__data.n_children);
          eval_data.r.fi = eval_data.r.f = 0.;
          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 0
      //**** THINK ABOUT THIS LATER !
      eval_data.r.f = __data.f_evaluate(-1, __data.params.nn(), *eval_data.z,
                                  *v_ind, eval_data.r.f, 0,
                                  &((*eval_data.islp_data)[__data.node_num]));
#endif
    }
  }

  void retrieve_from_cache(const expression_node& __data)
  { // cache is disabled
    // the parallel der_cache MUST be initialized, too
    if(__data.operator_type == EXPRINFO_LIN)
    {
      const matrix<double>::Row& lrw(eval_data.mod->lin[__data.params.nn()]);
      matrix<double>::Row::const_iterator _x, _e;

      eval_data.r.f = linalg_dot(lrw,*eval_data.z,0.);
      eval_data.r.fi = interval(0.);
      for(_x = lrw.begin(); _x != _e; ++_x)
        eval_data.r.fi += (*eval_data.z)[_x.index()] * interval(*_x);
    }
    else if(__data.operator_type == EXPRINFO_QUAD)
    {
      std::cerr << "Slope Evaluation of QUAD: NYI!" << std::endl;
      throw "NYI";
    }
    else
      eval_data.r.f = (*eval_data.f)[__data.node_num];
    eval_data.r.rg = (*eval_data.range)[__data.node_num];
  }

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

  int update(const expression_node& __data,
              const func_islp_return_type& __rval)
  {
    int ret = 0;
    interval __x;
    if(__data.operator_type < 0)
    {
      switch(__data.operator_type)
      {
        case EXPRINFO_CONSTANT:
          eval_data.r.rg = eval_data.r.f = __data.params.nd();
          // don't care about islp_data, CONSTANTS are ALWAYS leafs
          break;
        case EXPRINFO_VARIABLE:
          eval_data.r.fi = eval_data.r.f = (*eval_data.z)[__data.params.nn()];
          eval_data.r.rg = (*eval_data.range)[__data.node_num];
          // don't care about islp_data, VARIABLES are ALWAYS leafs
          break;
        case EXPRINFO_SUM:
        case EXPRINFO_MEAN:
          { double h = __data.coeffs[eval_data.n];
            eval_data.r.f += h*__rval.f;
            (*eval_data.islp_data)[__data.node_num][eval_data.n++] = h;
            eval_data.r.fi += h*(__rval.fi-__rval.f+interval(__rval.f));
          }
          break;
        case EXPRINFO_PROD:
          {
            (*eval_data.islp_data)[__data.node_num][eval_data.n] =
                                                                eval_data.r.rg;

            eval_data.r.fi *= __rval.f;
            eval_data.r.fi += eval_data.r.rg*(__rval.fi-__rval.f);

            eval_data.r.f *= __rval.f;
            eval_data.r.rg *= __rval.rg;

            for(int i = eval_data.n-1; i >= 0; i--)
              (*eval_data.islp_data)[__data.node_num][i] *= __rval.f;
          }
          ++eval_data.n;
          break;
        case EXPRINFO_MONOME:
          std::cerr << "func_islp_evaluator: MONOME NYI" << std::endl;
          throw "NYI";
#if 0
          if(eval_data.n == 0)
          {
            int n = __data.params.n()[0];
            if(n != 0)
            {
              eval_data.r = __power(__data.coeffs[0], __rval, n);
              (*eval_data.id_data)[__data.node_num][0] =
                    n*__power(__data.coeffs[0], __rval, n-1)*__data.coeffs[0];
            }
            else
            {
              eval_data.r = 1;
              (*eval_data.id_data)[__data.node_num][0] = 0;
            }
          }
          else
          {
            int n = __data.params.n()[eval_data.n];
            if(n != 0)
            {
              __x = __power(__data.coeffs[eval_data.n], __rval, n-1)*
                                __data.coeffs[eval_data.n];
              (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                                        eval_data.r*n*__x;
              __x = __power(__data.coeffs[eval_data.n], __rval, n);
              eval_data.r *= __x;
              for(int i = eval_data.n-1; i >= 0; i--)
                (*eval_data.id_data)[__data.node_num][i] *= __x;
            }
            else
              (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
          }
          ++eval_data.n;
          // intersect with the known range
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_MAX:
          // this is probably wrong (roundoff?) -> check this
          { double h = __rval.f * __data.coeffs[eval_data.n];
            if(h >= eval_data.r.f)
            {
              if(h > eval_data.r.f)
                (*eval_data.islp_data)[__data.node_num][eval_data.u.info] = 0;
              (*eval_data.islp_data)[__data.node_num][eval_data.n] = 
                                  __data.coeffs[eval_data.n];
              eval_data.r.f = h;
              eval_data.r.fi = h;
              eval_data.u.info = eval_data.n;
            }
            else
            {
              (*eval_data.islp_data)[__data.node_num][eval_data.n] = 0;
            }
          }
          ++eval_data.n;
          break;
        case EXPRINFO_MIN:
          // this is probably wrong (roundoff?) -> check this
          { double h = __rval.f * __data.coeffs[eval_data.n];
            if(h <= eval_data.r.f)
            {
              if(h < eval_data.r.f)
                (*eval_data.islp_data)[__data.node_num][eval_data.u.info] = 0;
              (*eval_data.islp_data)[__data.node_num][eval_data.n] =
                                  __data.coeffs[eval_data.n];
              eval_data.r.f = h;
              eval_data.r.fi = h;
              eval_data.u.info = eval_data.n;
            }
            else
            {
              (*eval_data.islp_data)[__data.node_num][eval_data.n] = 0;
            }
          }
          ++eval_data.n;
          break;
        case EXPRINFO_SCPROD:
          std::cerr << "func_islp_evaluator: SCPROD NYI" << std::endl;
          throw "NYI";
#if 0
          { interval h = __data.coeffs[eval_data.n]*__rval;
            // if the number of children is odd, the list is implicitly
            // lengthened by one null element (i.e. implicitly shortened)
            if(eval_data.n & 1) // y_n
            {
              eval_data.r += eval_data.u.d*h;
              (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                     eval_data.u.d*__data.coeffs[eval_data.n-1];
              (*eval_data.id_data)[__data.node_num][eval_data.n-1] =
                                                  h*__data.coeffs[eval_data.n];
            }
            else // x_n
              eval_data.u.d = h;
          }
          eval_data.n++;
          // intersect with the known range
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_NORM:
          std::cerr << "func_islp_evaluator: NORM NYI" << std::endl;
          throw "NYI";
#if 0
          { interval h = __data.coeffs[eval_data.n]*__rval;
            eval_data.r += sqr(h);
            (*eval_data.id_data)[__data.node_num][eval_data.n-1] =
                                                2*h*__data.coeffs[eval_data.n];
          }
          eval_data.n++;
          // intersect with the known range
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_INVERT:
          { interval __h;
            eval_data.r.f /= __rval.f;
            (*eval_data.islp_data)[__data.node_num][0] = __h =
                                                eval_data.r.f/__rval.rg;
            eval_data.r.fi /= __rval.f;
            eval_data.r.fi += __h*(__rval.fi-__rval.f);
          }
          break;
        case EXPRINFO_DIV:  // has two children
          if(eval_data.n++ == 0)
          {
            eval_data.r.f = __rval.f;
            eval_data.r.fi = __rval.fi-__rval.f;
          }
          else
          {
            double h = 1/__rval.f;
            interval __h, __k;
            __k = eval_data.r.f;
            eval_data.r.f *= __data.params.nd()*h;
            __h = __data.params.nd();
            __h /= __rval.f;
            (*eval_data.islp_data)[__data.node_num][0] = __h;
            eval_data.r.fi *= __h;
            __h *= -(*eval_data.range)[__data.node_num];
            (*eval_data.islp_data)[__data.node_num][1] = __h;
            eval_data.r.fi += __h*(__rval.fi-__rval.f);
            eval_data.r.fi += __data.params.nd()*__k/__rval.f;
          }
          break;
        case EXPRINFO_SQUARE:
          { interval __h;
            double h = __data.coeffs[0]*__rval.f+__data.params.nd();
            eval_data.r.f = h*h;
            (*eval_data.islp_data)[__data.node_num][0] = __h =
                           (h+__data.coeffs[0]*(__rval.rg+__data.params.nd()))*
                           __data.coeffs[0];
            eval_data.r.fi = __rval.f;
            eval_data.r.fi *= __data.coeffs[0];
            eval_data.r.fi += __data.params.nd();
            eval_data.r.fi = sqr(eval_data.r.fi);
            eval_data.r.fi += __h*(__rval.fi-__rval.f);
          }
          break;
        case EXPRINFO_INTPOWER:
          { int hl = __data.params.nn();
            if(hl == 0)
            {
              eval_data.r.f = 1;
              (*eval_data.islp_data)[__data.node_num][0] = 0;
              eval_data.r.fi = 1;
            }
            else
            {
              double kl = __data.coeffs[0]*__rval.f;
              interval km = __data.coeffs[0]*__rval.rg;
              interval _h_k = interval(__rval.f)*__data.coeffs[0];
              interval _h_h;
              switch(hl)
              {
                case 1:
                  eval_data.r.f = kl;
                  (*eval_data.islp_data)[__data.node_num][0] = __data.coeffs[0];
                  _h_h = 0;
                  break;
                case 2:
                  eval_data.r.f = kl*kl;
                  _h_h = (*eval_data.islp_data)[__data.node_num][0] =
                        __data.coeffs[0]*(kl+km);
                  _h_k = sqr(_h_k);
                  break;
                case -1:
                  eval_data.r.f = 1/kl;
                  _h_h = (*eval_data.islp_data)[__data.node_num][0] =
                                -__data.coeffs[0]/(kl*km);
                  _h_k = 1./_h_k;
                  break;
                case -2:
                  { // slope: range(-(x_+z)/(x_^2 z^2),-(x^+z)/(x^^2 z^2),
                    //              1/4z^3)
                    //        the last term if -2z in [x_,x^]
                    
                    double h = 1/kl;
                    double k = h*h;
                    eval_data.r.f = k;
                    interval _h;

                    _h = -((km.inf()+kl)*k)/(km.inf()*km.inf());
                    _h.hull(-((km.sup()+kl)*k)/(km.sup()*km.sup()));
                    if(_h.contains(-2*eval_data.r.f))
                      _h.hull(0.25/(h*k));
                    _h_h = (*eval_data.islp_data)[__data.node_num][0] =
                                                _h*__data.coeffs[0];
                    _h_k = 1./sqr(_h_k);
                  }
                  break;
                default:
                  if(hl > 0)
                  {
                    if(hl & 1) // odd
                    {
                      //****SUBOPTIMAL****
                      interval _h = (hl+0.0)*ipow(km,hl-1);
                      _h_h = (*eval_data.islp_data)[__data.node_num][0] =
                                                _h*__data.coeffs[0];
                    }
                    else // even
                    { // the function is convex
                      interval i1, i2;
                      double up;

                      i1 = km.sup();
                      if(kl != km.sup())
                      {
                        i1.ipow(hl);
                        i2 = kl;
                        i2.ipow(hl);
                        i1 -= i2;
                        i1 /= (km.sup()-kl);
                      }
                      else
                      {
                        i1.ipow(hl-1);
                        i1 *= hl;
                      }
                      up = i1.sup();

                      i1 = km.inf();
                      if(kl != km.inf())
                      {
                        i1.ipow(hl);
                        if(kl == km.sup())
                        { // otherwise i2 still contains kl^hl;
                          i2 = kl;
                          i2.ipow(hl);
                        }
                        i1 -= i2;
                        i1 /= (km.inf()-kl);
                      }
                      else
                      {
                        i1.ipow(hl-1);
                        i1 *= hl;
                      }
                      i1 = i1.inf();

                      _h_h = (*eval_data.islp_data)[__data.node_num][0] =
                                       __data.coeffs[0]*interval(i1.inf(),up);
                    }
                  }
                  else // hl < 0
                  {
                    if((-hl) & 1) // odd
                    {
                      //****SUBOPTIMAL****
                      interval _h = (hl+0.0)*ipow(km,hl-1);
                      _h_h = (*eval_data.islp_data)[__data.node_num][0] =
                                                _h*__data.coeffs[0];
                    }
                    else // even
                    {
                      //****SUBOPTIMAL****
                      interval _h = (hl+0.0)*ipow(km,hl-1);
                      _h_h = (*eval_data.islp_data)[__data.node_num][0] =
                                                _h*__data.coeffs[0];
                    }
                  }
                  _h_k = ipow(_h_k, hl);
                  break;
              }
              eval_data.r.fi = _h_h*(__rval.fi-__rval.f)+_h_k;
            }
          }
          break;
        case EXPRINFO_SQROOT:
          { double h = sqrt(__data.coeffs[0]*__rval.f+__data.params.nd());
            interval _h = sqrt(__data.coeffs[0]*__rval.rg+__data.params.nd());
            eval_data.r.f = h;
            _h = (*eval_data.islp_data)[__data.node_num][0] = __data.coeffs[0]/
              (h+_h);
            eval_data.r.fi = (__rval.fi-__rval.f)*_h+
              sqrt(__data.coeffs[0]*interval(__rval.f)+__data.params.nd());
          }
          break;
        case EXPRINFO_ABS:
          { double h = __data.coeffs[0]*__rval.f+__data.params.nd();
            interval _h = __data.coeffs[0]*__rval.rg+__data.params.nd();
            eval_data.r.f = fabs(h);
            if(_h.inf() >= 0 && h >= 0)
              _h = __data.coeffs[0];
            else if(_h.sup() <= 0 && h <= 0)
              _h = -__data.coeffs[0];
            else if(h > 0)
            {
              interval _lo,_up;
              if(_h.sup() < 0)
                _up = (interval(_h.sup())+h)/(h-_h.sup());
              else
                _up = 1.;
              _lo = (interval(_h.inf())+h)/(h-_h.inf());
              _h = __data.coeffs[0]*interval(_lo.inf(), _up.sup());
            }
            else if(h < 0)
            {
              interval _lo,_up;
              if(_h.inf() > 0)
                _lo = (interval(_h.inf())+h)/(h-_h.inf());
              else
                _lo = -1.;
              _up = (interval(_h.sup())+h)/(h-_h.sup());
              _h = __data.coeffs[0]*interval(_lo.inf(), _up.sup());
            }
            else
              _h = interval(-__data.coeffs[0],__data.coeffs[0]);
            (*eval_data.islp_data)[__data.node_num][0] = _h;
            eval_data.r.fi = (__rval.fi-__rval.f)*_h+
              abs(__data.coeffs[0]*interval(__rval.f)+__data.params.nd());
          }
          break;
        case EXPRINFO_POW:
          //****SUBOPTIMAL****
          if(eval_data.n++ == 0)
          {
            eval_data.r.f = __rval.f;
            eval_data.r.fi = __rval.fi-__rval.f;
            eval_data.r.rg = __rval.rg;
          }   
          else
          {
            double hh = __rval.f*__data.coeffs[1];
            if(hh == 0)
            {
              (*eval_data.islp_data)[__data.node_num][0] = 0;
              (*eval_data.islp_data)[__data.node_num][1] = 0;
              eval_data.r.f = 1;
              eval_data.r.fi = 1;
            }
            else
            {
              interval _hh = __rval.rg*__data.coeffs[0];
              interval _hf = interval(eval_data.r.f)*__data.coeffs[0]+
                                        __data.params.nd();

              eval_data.r.f *= __data.coeffs[0];
              eval_data.r.f += __data.params.nd();
              eval_data.r.rg *= __data.coeffs[0];
              eval_data.r.rg += __data.params.nd();

              double h = pow(eval_data.r.f, hh);
              interval _h = __rval.rg*pow(eval_data.r.rg, __rval.rg-1.)*
                              __data.coeffs[0];
              eval_data.r.fi *= _h;
              eval_data.r.fi += pow(_hf,interval(__rval.f)*__data.coeffs[1]);

              (*eval_data.islp_data)[__data.node_num][0] = _h;

              _h = (*eval_data.islp_data)[__data.node_num][1] =
                             log(eval_data.r.rg)*pow(eval_data.r.rg,__rval.rg)*
                             __data.coeffs[1];
              eval_data.r.f = h;
              eval_data.r.fi += (__rval.fi-__rval.f)*_h;
            }
          }
          break;
        case EXPRINFO_EXP:
          { double k = __rval.f*__data.coeffs[0]+__data.params.nd();
            double h = exp(k);
            interval _h = __rval.rg*__data.coeffs[0]+__data.params.nd();
            interval _lo(_h.inf());
            interval _up(_h.sup());
            _lo = (exp(_lo)-h)/(_lo-k);
            _up = (exp(_up)-h)/(_up-k);
            eval_data.r.f = h;
            _h = (*eval_data.islp_data)[__data.node_num][0] =
                interval(_lo.inf(),_up.sup())*__data.coeffs[0];
            eval_data.r.fi = (__rval.fi-__rval.f)*_h+
              exp(__data.coeffs[0]*interval(__rval.f)+__data.params.nd());
          }
          break;
        case EXPRINFO_LOG:
          { double k = __rval.f*__data.coeffs[0]+__data.params.nd();
            double h = log(k);
            interval _h = __rval.rg*__data.coeffs[0]+__data.params.nd();
            interval _lo(_h.sup());
            interval _up(_h.inf());
            _lo = (log(_lo)-h)/(_lo-k);
            _up = (log(_up)-h)/(_up-k);
            eval_data.r.f = h;
            _h = (*eval_data.islp_data)[__data.node_num][0] =
                interval(_lo.inf(),_up.sup())*__data.coeffs[0];
            eval_data.r.fi = (__rval.fi-__rval.f)*_h+
              log(__data.coeffs[0]*interval(__rval.f)+__data.params.nd());
          }
          break;
        case EXPRINFO_SIN:
          //****SUBOPTIMAL****
          { double k = __rval.f*__data.coeffs[0]+__data.params.nd();
            double h = sin(k);
            interval _h = __rval.rg*__data.coeffs[0]+__data.params.nd();
            eval_data.r.f = h;
            _h = (*eval_data.islp_data)[__data.node_num][0] =
                cos(_h)*__data.coeffs[0];
            eval_data.r.fi = (__rval.fi-__rval.f)*_h+
              sin(__data.coeffs[0]*interval(__rval.f)+__data.params.nd());
          }
          break;
        case EXPRINFO_COS:
          //****SUBOPTIMAL****
          { double k = __rval.f*__data.coeffs[0]+__data.params.nd();
            double h = cos(k);
            interval _h = __rval.rg*__data.coeffs[0]+__data.params.nd();
            eval_data.r.f = h;
            _h = (*eval_data.islp_data)[__data.node_num][0] =
                -sin(_h)*__data.coeffs[0];
            eval_data.r.fi = (__rval.fi-__rval.f)*_h+
              cos(__data.coeffs[0]*interval(__rval.f)+__data.params.nd());
          }
          break;
        case EXPRINFO_ATAN2:
          //****SUBOPTIMAL****
          if(eval_data.n++ == 0)
          {
            eval_data.r.f = __rval.f;
            eval_data.r.fi = __rval.fi-__rval.f;
            eval_data.r.rg = __rval.rg;
          }
          else
          {
            double hh = __rval.f * __data.coeffs[1];
            interval _hh = __rval.rg * __data.coeffs[1];
            eval_data.r.f *= __data.coeffs[0];
            eval_data.r.rg *= __data.coeffs[0]; 

            interval _h = sqr(eval_data.r.rg)+sqr(__rval.rg);
            interval _j;

            _j = (*eval_data.islp_data)[__data.node_num][0] =
                                                __data.coeffs[0]*_hh/_h;
            eval_data.r.fi *= _j;

            _j = (*eval_data.islp_data)[__data.node_num][1] =
                                           -__data.coeffs[1]*eval_data.r.rg/_h;
            eval_data.r.fi += _j*(__rval.fi*__data.coeffs[1]-hh);

            eval_data.r.fi += atan2(interval(eval_data.r.f),interval(hh));
            eval_data.r.f = atan2(eval_data.r.f, hh);
          }
          break;
        case EXPRINFO_GAUSS:
          //****SUBOPTIMAL****
          { double k = (__rval.f*__data.coeffs[0]-__data.params.d()[0])/
                                __data.params.d()[1];
            double h = exp(-k*k);
            interval _h = (__rval.rg*__data.coeffs[0]-__data.params.d()[0])/
                                __data.params.d()[1];
            eval_data.r.f = h;
            _h = (*eval_data.islp_data)[__data.node_num][0] =
                -2.*_h*exp(-sqr(_h))*__data.coeffs[0]/__data.params.d()[1];
            interval _k = (interval(__rval.f)*__data.coeffs[0]-
                        __data.params.d()[0])/__data.params.d()[1];
            eval_data.r.fi = (__rval.fi-__rval.f)*_h+exp(-sqr(_k));
          }
          break;
        case EXPRINFO_POLY:
          std::cerr << "Polynomes NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_LIN:
        case EXPRINFO_QUAD:
          // we will never be here
          break;
        case EXPRINFO_IN:
          std::cerr << "func_islp_evaluator: IN NYI" << std::endl;
          throw "NYI";
#if 0
          {
            __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(i.disjoint(__x))     // it is not in
            {
              eval_data.r = -1.;
              ret = -1;
              std::vector<interval>& r((*eval_data.id_data)[__data.node_num]);
              r.erase(r.begin(),r.end());
              r.insert(r.end(),__data.n_children,interval(0.));
            }
            else
            {
              if(!i.superset(__x))  // in and out are possible
              {
                if(i.sup() == __x.inf() || i.inf() == __x.sup())
                {
                  // intersection is only on the border
                  eval_data.r = interval(-1.,0.);
                  if(i.sup() == __x.inf())
                    (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                interval(-INFINITY,0.);
                  else
                    (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                interval(0.,INFINITY);
                }
                else
                {
                  // intersection is big
                  if(eval_data.r.sup() > 0)
                    eval_data.r = interval(-1.,1.);

                  if(i.subset(__x))
                    (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                interval(-INFINITY,INFINITY);
                  else if(__x.contains(i.inf()))
                    (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                interval(0.,INFINITY);
                  else
                    (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                interval(-INFINITY,0.);
                }
              }
              else if(i.inf() == __x.inf() || i.sup() == __x.sup())
              // they are not disjoint, and i is a superset of __x
              // but a border is possible
              {
                if(eval_data.r.inf() == 1.)
                {
                  if(__x.is_thin())       // only border
                    eval_data.r = 0.;
                  else                    // interior and border possible
                    eval_data.r = interval(0.,1.);
                }
                else if(__x.is_thin() && eval_data.r.inf() == 0.)
                  eval_data.r = 0.;
                if(__x.inf() == __x.sup())
                  (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
                else if(i.inf() == __x.inf() && i.sup() == __x.sup())
                  (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                  interval(-INFINITY,INFINITY);
                else if(i.inf() == __x.inf())
                  (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                  interval(0.,INFINITY);
                else
                  (*eval_data.id_data)[__data.node_num][eval_data.n] = 
                                                  interval(-INFINITY,0.);
              }
              else
              // if __x is in the interior of i eval_data.r does not change
                (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
            }
          }
          eval_data.n++;
#endif
          break;
        case EXPRINFO_IF:
          std::cerr << "func_islp_evaluator: IF NYI" << std::endl;
          throw "NYI";
#if 0
          __x = __rval * __data.coeffs[eval_data.n];
          if(eval_data.n == 0) // this is the condition
          {
            const interval& i(__data.params.ni());
            if(i.superset(__x))      // always true
            {
              eval_data.u.info = 0;
              (*eval_data.id_data)[__data.node_num][0] = 0;
              (*eval_data.id_data)[__data.node_num][2] = 0;
            }
            else if(i.disjoint(__x)) // never true
            {
              eval_data.u.info = 1;
              (*eval_data.id_data)[__data.node_num][0] = 0;
              (*eval_data.id_data)[__data.node_num][1] = 0;
              ret = 1;               // skip if part
            }
            else                     // both is possible
              eval_data.u.info = 1;
          }
          else if(eval_data.n == 1)  // this is either the if or the else part
          {
            eval_data.r = __x;
            (*eval_data.id_data)[__data.node_num][1] = __data.coeffs[1];
            if(eval_data.u.info == 0)  // it was always true
              ret = -1;                // don't continue walking
          }
          else // this is the else part if both possibilities could be true
          {
            if(eval_data.u.info == 1)
              eval_data.r = __x;
            else
            {
              eval_data.r = eval_data.r.hull(__x);
              (*eval_data.id_data)[__data.node_num][0] =
                                                interval(-INFINITY,INFINITY);
            }
            (*eval_data.id_data)[__data.node_num][2] = __data.coeffs[2];
          }
          if(eval_data.n == 2 || ret == -1)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
          else
            eval_data.n += 1+ret;
#endif
          break;
        case EXPRINFO_AND:
          std::cerr << "func_islp_evaluator: AND NYI" << std::endl;
          throw "NYI";
#if 0
          { __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(i.disjoint(__x))
            {
              eval_data.r = 0;
              std::vector<interval>& r((*eval_data.id_data)[__data.node_num]);
              r.erase(r.begin(),r.end());
              r.insert(r.end(),__data.n_children,interval(0.));
              ret = -1;
            }
            else if(!i.superset(__x))
            {
              double lo = 0, hi = 0;
              if(eval_data.r != 0)
                eval_data.r = interval(0.,1.);
              if(__x.inf() < i.inf())
                hi = INFINITY;
              if(__x.sup() > i.sup())
                lo = -INFINITY;
              (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                                               interval(lo,hi);
            }
            else
              (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
          }
          eval_data.n++;
          // intersect with the known range
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_OR:
          std::cerr << "func_islp_evaluator: OR NYI" << std::endl;
          throw "NYI";
#if 0
          { __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(i.superset(__x))
            {
              eval_data.r = 1;
              std::vector<interval>& r((*eval_data.id_data)[__data.node_num]);
              r.erase(r.begin(),r.end());
              r.insert(r.end(),__data.n_children,interval(0.));
              ret = -1;
            }
            else if(!i.disjoint(__x))
            {
              if(eval_data.r != 0)
                eval_data.r = interval(0.,1.);
              double lo = 0, hi = 0;
              if(__x.inf() < i.inf())
                hi = INFINITY;
              if(__x.sup() > i.sup())
                lo = -INFINITY;
              (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                                               interval(lo,hi);
            }
            else
              (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
          }
          eval_data.n++;
          // intersect with the known range
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_NOT:
          std::cerr << "func_islp_evaluator: NOT NYI" << std::endl;
          throw "NYI";
#if 0
          { __x = __data.coeffs[0]*__rval;
            const interval& i(__data.params.ni());
            if(i.superset(__x))
            {
              eval_data.r = 0;
              (*eval_data.id_data)[__data.node_num][0] = 0;
            }
            else if(i.disjoint(__x))
            {
              eval_data.r = 1;
              (*eval_data.id_data)[__data.node_num][0] = 0;
            }
            else
            {
              eval_data.r = interval(0.,1.);
              double lo = 0, hi = 0;
              if(__x.inf() < i.inf())
                lo = -INFINITY;
              if(__x.sup() > i.sup())
                hi = INFINITY;
              (*eval_data.id_data)[__data.node_num][0] = interval(lo,hi);
            }
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
          }
#endif
          break;
        case EXPRINFO_IMPLIES:
          std::cerr << "func_islp_evaluator: IMPLIES NYI" << std::endl;
          throw "NYI";
#if 0
          { const interval& i(__data.params.i()[eval_data.n]);
            __x = __rval * __data.coeffs[eval_data.n];
            if(eval_data.n == 0)
            {
              if(i.disjoint(__x))
              {
                eval_data.r = 1.;
                (*eval_data.id_data)[__data.node_num][0] = 0;
                (*eval_data.id_data)[__data.node_num][1] = 0;
                ret = -1;        // ex falso quodlibet
              }
              else if(!i.superset(__x))
              {
                eval_data.r = interval(0.,1.);
                double lo = 0, hi = 0;
                if(__x.inf() < i.inf())
                  lo = -INFINITY;
                if(__x.sup() > i.sup())
                  hi = INFINITY;
                (*eval_data.id_data)[__data.node_num][0] = interval(lo,hi);
              }
              else
              {
                eval_data.r = interval(0.);
                (*eval_data.id_data)[__data.node_num][0] = 0;
              }
            }
            else // we only get here if first clause might be true
            {
              if(i.superset(__x)) // true result is always true
              {
                eval_data.r = 1.;
                (*eval_data.id_data)[__data.node_num][0] = 0;
                (*eval_data.id_data)[__data.node_num][1] = 0;
              }
              else if(!i.disjoint(__x))
              {
                eval_data.r = interval(0.,1.);
                double lo = 0, hi = 0;
                if(__x.inf() < i.inf())
                  hi = INFINITY;
                if(__x.sup() > i.sup())
                  lo = -INFINITY;
                (*eval_data.id_data)[__data.node_num][0] = interval(lo,hi);
              }
              else
                (*eval_data.id_data)[__data.node_num][1] = 0;
              eval_data.r.intersect((*eval_data.range)[__data.node_num]);
            }
            ++eval_data.n;
          }
#endif
          break;
        case EXPRINFO_COUNT:
          std::cerr << "func_islp_evaluator: COUNT NYI" << std::endl;
          throw "NYI";
#if 0
          { __x = __data.coeffs[eval_data.n]*__rval;
            const interval& i(__data.params.i()[eval_data.n]);
            if(i.superset(__x))
            {
              eval_data.r += 1;
              (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
            }
            else if(!i.disjoint(__x))
            {
              eval_data.r += interval(0.,1.);
              double lo = 0, hi = 0;
              if(__x.inf() < i.inf())
                hi = INFINITY;
              if(__x.sup() > i.sup())
                lo = -INFINITY;
              (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                                               interval(lo,hi);
            }
            else
              (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
          }
          eval_data.n++;
          // intersect with the known range
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_ALLDIFF:
          std::cerr << "func_islp_evaluator: ALLDIFF NYI" << std::endl;
          throw "NYI";
#if 0
          { int i; 
            std::vector<interval>::const_iterator _b;
            __x = __data.coeffs[eval_data.n]*__rval;
            (*eval_data.id_data)[__data.node_num][eval_data.n] = 0;
            for(_b = ((std::vector<interval>*)eval_data.u.p)->begin(), i = 0;
                _b != ((std::vector<interval>*)eval_data.u.p)->end(); ++_b, ++i)
            {
              interval __h = abs(__x-*_b);
              if(__h.sup() <= __data.params.nd())  // they are always equal
              {
                eval_data.r = 0;
                std::vector<interval>& r((*eval_data.id_data)[__data.node_num]);
                r.erase(r.begin(),r.end());
                r.insert(r.end(),__data.n_children,interval(0.));
                ret = -1;               // don't continue
                break;
              }
              else if(__h.inf() <= __data.params.nd()) // they are not always different
              {
                eval_data.r = interval(0.,1.);
                (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                                interval(-INFINITY,INFINITY);
                (*eval_data.id_data)[__data.node_num][i] =
                                                interval(-INFINITY,INFINITY);
              }
            }
            if(ret != -1)
              ((std::vector<interval>*) eval_data.u.p)->push_back(__x);
          }
          eval_data.n++;
          if(eval_data.n == __data.n_children || ret == -1)
          {
            if(eval_data.r == 1.)       // all are always different
            {
              std::vector<interval>& r((*eval_data.id_data)[__data.node_num]);
              r.erase(r.begin(),r.end());
              r.insert(r.end(),__data.n_children,interval(0.));
            }
            delete (std::vector<interval>*) eval_data.u.p;
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
          }
#endif
          break;
        case EXPRINFO_HISTOGRAM:
          std::cerr << "func_islp_evaluator: histogram NYI" << std::endl;
          throw "NYI";
        case EXPRINFO_LEVEL:
          std::cerr << "func_islp_evaluator: LEVEL NYI" << std::endl;
          throw "NYI";
#if 0
          { // we have to make this one more complicated than
            // in the interval evaluation. For deciding the final
            // outcome, we (seemingly) need to know the results
            // for all before determining the derivatives
            int lo = 0;
            int hi = 0;
            __x = __data.coeffs[eval_data.n]*__rval;
            interval _h;

            // first adjust upper bound
            if(hi != INT_MAX)
            {
              while(hi < __data.params.im().nrows())
              { 
                _h = __data.params.im()[hi][eval_data.n];
                if(_h.superset(__x))
                  break;
                hi++;
              }
              if(hi == __data.params.im().nrows())
                hi = INT_MAX;
            }
            // then adjust lower bound
            if(lo != INT_MAX)
            {
              while(lo < __data.params.im().nrows())
              { 
                _h = __data.params.im()[lo][eval_data.n];
                if(!_h.disjoint(__x))
                  break;
                lo++;
              }
              if(lo == __data.params.im().nrows())
                lo = INT_MAX;
            }
            ((std::vector<interval>*)eval_data.u.p)->
                                        push_back(interval(lo+0.,hi+0.));
            eval_data.r = interval(min(eval_data.r.inf(),lo+0.),
                                   max(eval_data.r.sup(),hi+0.));
          }
          eval_data.n++;
          if(eval_data.n == __data.n_children)
          {
            std::vector<interval>& v(*((std::vector<interval>*)eval_data.u.p));
            for(int i = 0; i < __data.n_children; ++i)
            {
              if(v.is_thin() || v.sup() < eval_data.r.inf())
                (*eval_data.id_data)[__data.node_num][i] = 0;
              else
                (*eval_data.id_data)[__data.node_num][i] =
                                                interval(0.,INFINITY);
            }
            delete (std::vector<interval>*)eval_data.u.p;
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
          }
#endif
          break;
        case EXPRINFO_NEIGHBOR:
          std::cerr << "func_islp_evaluator: NEIGHBOR NYI" << std::endl;
          throw "NYI";
#if 0
          // this is suboptimal. It does not take into account
          // all the possibilities to proove that the result actually
          // EQUALS 1
          // furthermore, it sets all derivatives to [-I,I] because
          // computing derivatives is nonsense, anyway
          (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                                interval(-INFINITY,INFINITY);
          if(eval_data.n == 0)
            eval_data.r = __data.coeffs[0]*__rval;
          else
          {
            interval h = eval_data.r;
            eval_data.r = 0;
            __x = __data.coeffs[1]*__rval;
            for(unsigned int i = 0; i < __data.params.n().size(); i+=2)
            {
              if(h == __data.params.n()[i]+0. &&
                 __x == __data.params.n()[i+1]+0.)
              {
                eval_data.r = 1;
                break;
              }
              else if(h.contains(__data.params.n()[i]+0.) &&
                      __x.contains(__data.params.n()[i+1]+0.))
              {
                eval_data.r = interval(0.,1.);
                break;
              }
            }
          }
          eval_data.n++;
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_NOGOOD:
          std::cerr << "func_islp_evaluator: NOGOOD NYI" << std::endl;
          throw "NYI";
#if 0
          // derivatives are nonsense here, so set them to [-I,I]
          (*eval_data.id_data)[__data.node_num][eval_data.n] =
                                                interval(-INFINITY,INFINITY);
          __x = __data.coeffs[eval_data.n]*__rval;
          if(eval_data.r != 0.)
          { 
            if(!__x.contains(__data.params.n()[eval_data.n]+0.))
            {
              eval_data.r = 0.;
              ret = -1;
            }
            else if(__x != __data.params.n()[eval_data.n]+0.)
              eval_data.r = interval(0.,1.);
          }
          eval_data.n++;
          if(eval_data.n == __data.n_children)
            eval_data.r.intersect((*eval_data.range)[__data.node_num]);
#endif
          break;
        case EXPRINFO_EXPECTATION:
          std::cerr << "func_islp_evaluator: E NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_INTEGRAL:
          std::cerr << "func_islp_evaluator: INT NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_LOOKUP:
        case EXPRINFO_PWLIN:
        case EXPRINFO_SPLINE:
        case EXPRINFO_PWCONSTLC:
        case EXPRINFO_PWCONSTRC:
          std::cerr << "func_islp_evaluator: Table Operations NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_DET:
        case EXPRINFO_COND:
        case EXPRINFO_PSD:
        case EXPRINFO_MPROD:
        case EXPRINFO_FEM:
          std::cerr << "func_islp_evaluator: Matrix Operations NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_RE:
        case EXPRINFO_IM:
        case EXPRINFO_ARG:
        case EXPRINFO_CPLXCONJ:
          std::cerr << "func_islp_evaluator: Complex Numbers NYI" << std::endl;
          throw "NYI";
          break;
        case EXPRINFO_CMPROD:
        case EXPRINFO_CGFEM:
          std::cerr << "func_islp_evaluator: Const Matrix Operations NYI" << std::endl;
          throw "NYI";
          break;
        default:
          std::cerr << "func_islp_evaluator: Unknown operator_type "
                << __data.operator_type << "!" << std::endl;
          throw "Programming Error";
          // give bad messages
          break;
      }
    }
    else if(__data.operator_type > 0)
    {
#if 0
      //***** THINK HERE LATER
      // update the function arguments
      eval_data.r.f = __data.f_evaluate(eval_data.n++, __data.params.nn(),
                        *eval_data.z, *v_ind, eval_data.r.f, __rval,
                            &(*eval_data.islp_data)[__data.node_num]);
#endif
    }
    // keep the centers
    (*eval_data.f)[__data.node_num] = eval_data.r.f;
    if(__data.n_children == 1 || (eval_data.n == __data.n_children &&
          __data.n_children > 0))
      eval_data.r.rg = (*eval_data.range)[__data.node_num];
    return ret;
  }

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

struct islp_eval_type
{
  std::vector<std::vector<interval> >* islp_data;
  std::vector<std::vector<interval> >* islp_cache;
  std::vector<interval>* islp_vec;
  const model* mod;
  interval mult;
  interval mult_trans;
  unsigned int child_n;
};

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

protected:
  bool is_cached(const node_data_type& __data)
     { 
       if(eval_data.islp_cache && __data.n_parents > 1 && __data.n_children > 0
          && (*eval_data.islp_cache)[__data.node_num].size() > 0 &&
          v_ind->match(__data.var_indicator()))
       {
         return true;
       }
       else
         return false;
     }
                                        
public:
  islp_eval(std::vector<std::vector<interval> >& __der_data,
            variable_indicator& __v, const model& __m,
            std::vector<std::vector<interval > >* __d, 
            std::vector<interval>& __islp)
  { // __islp MUST BE INITIALIZED WITH 0 AND THE CORRECT LENGTH
    eval_data.islp_data = &__der_data;
    eval_data.islp_cache = __d;
    eval_data.islp_vec = &__islp;
    eval_data.mod = &__m;
    eval_data.mult_trans = 1;
    eval_data.mult = 0;
    v_ind = &__v;
  }

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

  ~islp_eval() {}

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

  void new_result(std::vector<interval>& __islp)
  {
    eval_data.islp_vec = &__islp;
  }

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

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

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

  // before the children are visited
  int calculate(const expression_node& __data)
  {
    if(__data.operator_type == EXPRINFO_CONSTANT)
      return 0;
    else if(__data.operator_type == EXPRINFO_VARIABLE)
    {
      // update the slope
      (*eval_data.islp_vec)[__data.params.nn()] += eval_data.mult_trans;
      return 0;
    }
    else if(__data.operator_type == EXPRINFO_LIN)
    {
      const matrix<double>::Row& lrw(eval_data.mod->lin[__data.params.nn()]);
      matrix<double>::Row::const_iterator _x, _e;
      _e = lrw.end();
      for(_x = lrw.begin(); _x != _e; ++_x)
        (*eval_data.islp_vec)[_x.index()] += *_x * eval_data.mult_trans;
      return 0;     // don't perform the remaining graph walk
    }
    else if(__data.operator_type == EXPRINFO_QUAD)
    {
      std::cerr << "islp_evaluator: QUAD is NYI!" << std::endl;
      throw "NYI";
      return 0;     // don't perform the remaining graph walk
    }
    else if(__data.ev && (*__data.ev)[DER_EVALUATOR])
    // is there a short-cut evaluator defined?
    {
      std::vector<interval> __h = (*(islp_evaluator)(*__data.ev)[ISLP_EVALUATOR])(
                              (*eval_data.islp_data)[__data.node_num], *v_ind);

      for(unsigned int i = 0; i <= (*eval_data.islp_vec).size(); ++i)
        (*eval_data.islp_vec)[i] += eval_data.mult*__h[i];
      return 0;     // don't perform the remaining graph walk
    }
    else if(eval_data.mult_trans.inf() == 0 && eval_data.mult_trans.sup() == 0)
    // does not contribute to islp
      return 0;
    else
    {
      eval_data.child_n = 1;
      eval_data.mult = eval_data.mult_trans;
      if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.islp_cache)
      { 
        eval_data.mult_trans = (*eval_data.islp_data)[__data.node_num][0];
      }
      else // set multiplier for first child
        eval_data.mult_trans *= (*eval_data.islp_data)[__data.node_num][0];
      return 1;
    }
  }

  // after all children have finished
  void cleanup(const expression_node& __data)
  {
    // use the cache only if it is worthwhile
    if(__data.n_parents > 1 && __data.n_children > 0 && eval_data.islp_cache
       && (*eval_data.islp_cache)[__data.node_num].size() == 0)
    {
      (*eval_data.islp_cache)[__data.node_num] = *eval_data.islp_vec;
      for(unsigned int i = 0; i <= (*eval_data.islp_vec).size(); ++i)
        (*eval_data.islp_vec)[i] *= eval_data.mult;
    }
  }

  void retrieve_from_cache(const expression_node& __data)
  {
    // the parallel f_cache MUST be initialized, too
    for(unsigned int i = 0; i <= (*eval_data.islp_vec).size(); ++i)
      (*eval_data.islp_vec)[i] = eval_data.mult_trans *
          (*eval_data.islp_cache)[__data.node_num][i];
  }

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

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

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

#endif /* _ISLP_EVALUATOR_H_ */

---