---

model-inline.h

Go to the documentation of this file.

// Model 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 _MODEL_INLINE_H_
#define _MODEL_INLINE_H_

#define MODEL_INLINE_DEBUG 0
#define MODEL_INLINE_DEBUG_SIMPLIFY 0

inline bool model_iddata::delete_ref(model_gid& __m)
{
  for(std::list<model_gid*>::iterator __w = backref.begin(); __w != backref.end();
      ++__w)
  {
    if(*__w == &__m)
    {
      backref.erase(__w);
      break;
    }
  }
  return backref.empty();
}

inline unsigned int model_iddata::get_node_id()
{
  if(node_free.size() > 0)
  {
    unsigned int ___a = node_free.back();
    node_free.pop_back();
    return ___a;
  }
  else
    return node_num_max++;
}

inline void model_iddata::remove_node_id(unsigned int _n)
{
  if(_n == node_num_max-1)
    --node_num_max;
  else
    node_free.push_back(_n);
  for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
      ++__b)
    (*__b)->remove_node_ref(_n);
}

inline unsigned int model_iddata::get_var_id()
{
  if(var_free.size() > 0)
  {
    unsigned int ___a = var_free.back();
    var_free.pop_back();
    return ___a;
  }
  else
    return var_num_max++;
}

inline void model_iddata::remove_var_id(unsigned int _n)
{
  if(_n == var_num_max-1)
    --var_num_max;
  else
    var_free.push_back(_n);
  for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
      ++__b)
    (*__b)->remove_var_ref(_n);
  if(var_names.size() > _n)
    var_names[_n] = std::string();
}

inline unsigned int model_iddata::get_const_id()
{
  if(const_free.size() > 0)
  {
    unsigned int ___a = const_free.back();
    const_free.pop_back();
    return ___a;
  }
  else
    return const_num_max++;
}

inline void model_iddata::remove_const_id(unsigned int _n)
{
  if(_n == const_num_max-1)
    --const_num_max;
  else
    const_free.push_back(_n);
  for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
      ++__b)
    (*__b)->remove_const_ref(_n);
  if(const_names.size() > _n)
    const_names[_n] = std::string();
}

inline void model_iddata::number_of_nodes(unsigned int _n)
{
  node_num_max = _n;
  for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
      ++__b)
    (*__b)->upd_number_of_nodes(_n);
}

inline void model_iddata::number_of_variables(unsigned int _n)
{
  var_num_max = _n;
  var_names.insert(var_names.end(), _n - var_names.size(), std::string());
  for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
      ++__b)
    (*__b)->upd_number_of_variables(_n);
}

inline void model_iddata::number_of_constraints(unsigned int _n)
{
  const_num_max = _n;
  const_names.insert(const_names.end(), _n - const_names.size(), std::string());
  for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
      ++__b)
    (*__b)->upd_number_of_constraints(_n);
}

inline const std::string model_iddata::var_name(unsigned int n) const
{
  if(n >= var_names.size() || var_names[n] == std::string())
  {
    char __s[32];
    sprintf(__s, "%%x%d", n);
    return std::string(__s);
  }
  else
   return var_names[n];
}

inline void model_iddata::var_name(unsigned int n, const std::string& vn)
{
  if(n >= var_names.size())
  {
    if(n > var_names.size())
      var_names.insert(var_names.end(), n-var_names.size(), std::string());
    var_names.push_back(vn);
  }
  else
    var_names[n] = vn;
}

inline const std::string model_iddata::const_name(unsigned int n) const
{
  if(n >= const_names.size() || const_names[n] == std::string())
  {
    char __s[32];
    sprintf(__s, "%%c%d", n);
    return std::string(__s);
  }
  else
   return const_names[n];
}

inline void model_iddata::const_name(unsigned int n, const std::string& vn)
{
  if(n >= const_names.size())
  {
    if(n > const_names.size())
      const_names.insert(const_names.end(), n-const_names.size(), std::string());
    const_names.push_back(vn);
  }
  else
    const_names[n] = vn;
}

inline const std::string model_iddata::obj_name() const
{
  if(obj_names == std::string())
    return std::string("%o");
  else
    return obj_names;
}

inline void model_iddata::obj_name(const std::string& vn) { obj_names = vn; }
inline double model_iddata::obj_adj() const { return obj_adjs; }
inline void model_iddata::obj_adj(double adj) { obj_adjs = adj; }
inline double model_iddata::obj_mult() const { return obj_mults; }
inline void model_iddata::obj_mult(double mult) { obj_mults = mult; }

inline std::pair<const std::string, double> model_iddata::fixed_var( unsigned int n)
  const
{
  if(n >= fixed_vars.size())
  {
    std::cerr << "index error in fixed_var!" << std::endl;
    throw "Programming error";
  }
  return fixed_vars[n];
}

inline void model_iddata::fixed_var(const std::string& vn, double val)
{
  fixed_vars.push_back(std::make_pair(vn, val));
}

inline const std::string& model_iddata::unused_var(unsigned int n) const
{
  if(n >= unused_vars.size())
  {
    std::cerr << "index error in unused_var!" << std::endl;
    throw "Programming error";
  }
  return unused_vars[n];
}

inline void model_iddata::unused_var(const std::string& vn)
{
  unused_vars.push_back(vn);
}

inline const std::string& model_iddata::unused_constr(unsigned int n) const
{
  if(n >= unused_constrs.size())
  {
    std::cerr << "index error in unused_constr!" << std::endl;
    throw "Programming error";
  }
  return unused_constrs[n];
}

inline void model_iddata::unused_constr(const std::string& vn)
{
  unused_constrs.push_back(vn);
}

inline void model_gid::upd_number_of_nodes(unsigned int _n)
{
  int __x;
  if((__x = _n - glob_ref.size()) > 0)
    glob_ref.insert(glob_ref.end(), __x, model_ref.ground());
}

inline void model_gid::upd_number_of_variables(unsigned int _n)
{
  int __x;
  unsigned int _n_old = gvar_ref.size();

  if((__x = _n - gvar_ref.size()) > 0)
  {
    model::walker _gr = model_ref.ground();

    gvar_ref.insert(gvar_ref.end(), __x, _gr);
    for(model::ref_iterator __b = glob_ref.begin(); __b != glob_ref.end();
        ++__b)
    {
      if((*__b) != _gr)
      {
        (*__b)->v_i.reserve(_n);
        (*__b)->v_i.unset(_n_old, _n-1);
      }
    }
  }
}

inline void model_gid::upd_number_of_constraints(unsigned int _n)
{
  int __x;
  if((__x = _n - gconst_ref.size()) > 0)
    gconst_ref.insert(gconst_ref.end(), __x, model_ref.ground());
}

inline void model_gid::renumber_variables(
    std::vector<std::pair<unsigned int,unsigned int> >& __v)
{
  int i;
  for(std::vector<std::pair<unsigned int,unsigned int> >::iterator __i = __v.begin();
      __i != __v.end(); ++__i)
  {
    i = (*__i).second;
    gvar_ref[i] = gvar_ref[(*__i).first];
    gvar_ref[i]->params = i;
  }
}

inline void model_gid::renumber_constraints(
    std::vector<std::pair<unsigned int,unsigned int> >& __v)
{
  unsigned int i;
  std::map<unsigned int, unsigned int>::iterator __mf;
  for(std::vector<std::pair<unsigned int,unsigned int> >::iterator __i =
      __v.begin(); __i != __v.end(); ++__i)
  {
    i = (*__i).second;
    gconst_ref[i] = gconst_ref[(*__i).first];
    unsigned int n = gconst_ref[i]->node_num;
    __mf = const_back_ref.find(n);
    if(__mf == const_back_ref.end())
      throw "This can't happen in renumber_constraints!";
    const_back_ref.erase(__mf);
    const_back_ref.insert(std::make_pair(n, i));
  }
}

inline void model_gid::mk_globref(unsigned int n, const model::walker& __w)
{
  if(glob_ref.size() > n)
    glob_ref[n] = __w;
  else if(n == glob_ref.size())
    glob_ref.push_back(__w);
  else
  {
    glob_ref.insert(glob_ref.end(), n-glob_ref.size()+1, model_ref.ground());
    glob_ref[n] = __w;
  }
}

inline void model_gid::mk_gvarref(unsigned int n, const model::walker& __w)
{
  if(gvar_ref.size() > n)
    gvar_ref[n] = __w;
  else if(n == gvar_ref.size())
    gvar_ref.push_back(__w);
  else
  {
    gvar_ref.insert(gvar_ref.end(), n-gvar_ref.size()+1, model_ref.ground());
    gvar_ref[n] = __w;
  }
}

inline void model_gid::mk_gconstref(unsigned int n, const model::walker& __w)
{
  if(n == gconst_ref.size())
    gconst_ref.push_back(__w);
  else
  {
    if(gconst_ref.size() <= n)
      gconst_ref.insert(gconst_ref.end(), n-gconst_ref.size()+1,
                      model_ref.ground());
    gconst_ref[n] = __w;
    std::map<unsigned int, unsigned int>::iterator __mf =
      const_back_ref.find(__w->node_num);
    if(__mf != const_back_ref.end())
      const_back_ref.erase(__mf);
  }
  const_back_ref.insert(std::make_pair(__w->node_num, n));
}

inline void model_gid::remove_const_ref(unsigned int _n)
{
  std::map<unsigned int, unsigned int>::iterator __mf;
  if(gconst_ref.size() > _n)
  {
    __mf = const_back_ref.find(gconst_ref[_n]->node_num);
    if(__mf != const_back_ref.end())
      const_back_ref.erase(__mf);
    if(gconst_ref.size() == _n) gconst_ref.pop_back();
    else gconst_ref[_n] = model_ref.ground();
  }
  else
  {
    for(__mf = const_back_ref.begin(); __mf != const_back_ref.end(); ++__mf)
      if((*__mf).second == _n)
      {
        const_back_ref.erase(__mf);
        break;
      }
  }
}

inline bool model_gid::get_const_num(unsigned int node_num,
                                     unsigned int& const_num)
{
  std::map<unsigned int, unsigned int>::iterator __mf =
                                        const_back_ref.find(node_num);
  if(__mf == const_back_ref.end())
    return false;
  const_num = (*__mf).second;
  return true;
}

inline void model_gid::make_const_back_ref(unsigned int node, unsigned int cnum)
{
  std::map<unsigned int, unsigned int>::iterator __mf =
                                        const_back_ref.find(node);
  if(__mf != const_back_ref.end())
    const_back_ref.erase(__mf);
  const_back_ref.insert(std::make_pair(node, cnum));
}

inline bool model::is_empty(const walker& __w) const { return gid->empty(__w); }
inline model::walker model::empty_reference() const
        { return gid->empty_reference(); }

inline void model::detach_gid() { gid = new model_gid(*this, *gid); }

inline model::model(int __num_of_vars) : _Base(), node_ref(), var_ref(),
                                ghost_ref(), keep_gid(true), lin(), matd(),
                                mati(), ocoeff(0), objective(ground()),
                                constraints()
{ gid = new model_gid(*this, __num_of_vars); }

inline model::model(model_gid* __id, bool clone) :
                                       _Base(), node_ref(), var_ref(),
                                       ghost_ref(), lin(), matd(), mati(),
                                       ocoeff(0), objective(ground()),
                                       constraints()
{
  keep_gid = true;
  if(__id) 
  {
    if(clone)
    {
      gid = __id;
      keep_gid = false;
    }
    else
      gid = new model_gid(*this, *__id);
  }
  else
    gid = new model_gid(*this);
}

inline model::model(const model& __m) : _Base(__m), node_ref(),
                                 var_ref(), ghost_ref(), gid(), keep_gid(true),
                                 lin(__m.lin), matd(__m.matd),
                                 mati(__m.mati), ocoeff(__m.ocoeff),
                                 constraints()
{
#if MODEL_INLINE_DEBUG
  std::cout << "Copy constructor of model called " << std::endl;
#endif
  gid = new model_gid(*this, *__m.gid_data());

  copy_contents(gid, __m);
}

inline model::model(model_gid* __id, const model& __m) : _Base(__m),
                                 node_ref(), var_ref(), ghost_ref(), gid(__id),
                                 keep_gid(false),
                                 lin(__m.lin), matd(__m.matd),
                                 mati(__m.mati), ocoeff(__m.ocoeff),
                                 constraints()
{
  if(!__id)
  {
    std::cerr << "NULL gid is forbidden in model constructor!" <<  std::endl;
    throw "Programming error";
  } 
  copy_contents(gid, __m);
}

inline model::~model()
{
#if MODEL_INLINE_DEBUG
  // brice 10/10/2002 : added a warning
  std::cerr << "Destructor of model called" << std::endl;
  if(keep_gid)
    delete gid;
#endif
}

inline bool model::get_const_num(unsigned int node_num,
                                 unsigned int& const_num) const
        { return gid->get_const_num(node_num, const_num); }
 

inline const std::string model::var_name(unsigned int n) const
        { return gid->var_name(n); }
inline const std::string model::const_name(unsigned int n) const
        { return gid->const_name(n); }
inline const std::string model::obj_name() const { return gid->obj_name(); }
inline double model::obj_adj() const { return gid->obj_adj(); }
inline double model::obj_mult() const { return gid->obj_mult(); }
inline size_t model::n_fixed_vars() const { return gid->n_fixed_vars(); }
inline std::pair<const std::string, double> model::fixed_var(unsigned int n) const
        { return gid->fixed_var(n); }
inline size_t model::n_unused_vars() const { return gid->n_unused_vars(); }
inline const std::string& model::unused_var(unsigned int n) const
        { return gid->unused_var(n); }
inline size_t model::n_unused_constrs() const
        { return gid->n_unused_constrs(); }
inline const std::string& model::unused_constr(unsigned int n) const
        { return gid->unused_constr(n); }
                                            
inline unsigned int model::number_of_variables() const
{ return gid->number_of_variables(); }

inline unsigned int model::number_of_constraints() const
{ return gid->number_of_constraints(); }

inline unsigned int model::number_of_nodes() const
{ return gid->number_of_nodes(); }

inline const model::walker& model::node(unsigned int i) const
{ return gid->node(i); }

inline const model::walker& model::var(unsigned int i) const
{ return gid->variable(i); }

inline const model::walker& model::constraint(unsigned int i) const
{ return gid->constraint(i); }

class model::sort_constraints :
        public std::binary_function<walker, walker, bool>
{
public:
  bool operator()(const walker& _c1, const walker& _c2) const
  {
    if(_c1.is_leaf() != _c2.is_leaf())
      if(_c1.is_leaf()) return true;
      else if(_c2.is_leaf()) return false;
    if(_c1->sem.degree != _c2->sem.degree)
      if(_c1->sem.degree == 1) return true;
      else if(_c2->sem.degree == 1) return false;
    if(_c1->f_bounds.is_thin() != _c2->f_bounds.is_thin())
      return _c1->f_bounds.is_thin() ? true : false;
    else if(_c1->f_bounds.is_bounded() == _c2->f_bounds.is_bounded() ||
      _c1->f_bounds.is_unbounded_below() == _c2->f_bounds.is_unbounded_below()||
      _c1->f_bounds.is_unbounded_above() == _c2->f_bounds.is_unbounded_above())
      return _c1->node_num < _c2->node_num;
    else if(_c1->f_bounds.is_unbounded_below())
      return true;
    else if(_c2->f_bounds.is_unbounded_below())
      return false;
    else 
      return _c1->f_bounds.is_bounded() ? false : true;
  }
};

#define SIMPLIFY_0_IS_CONST             1
#define SIMPLIFY_0_CONST_IS_INTEGER     (1<<24)

#define SIMPLIFY_0_IS_VAR               (1<<1)

#define SIMPLIFY_0_IS_SUM               (1<<2)
#define SIMPLIFY_0_SUM_IS_SIMPLE        (1<<24)

#define SIMPLIFY_0_IS_PROD              (1<<3)
#define SIMPLIFY_0_PROD_IS_SIMPLE       (1<<24)

#define SIMPLIFY_0_IS_MULTIPLICATIVE    (1<<4)
#define SIMPLIFY_0_IS_CORRECTEDMULT     (1<<5)

#define SIMPLIFY_0_IS_GHOST             (1<<31)

class model::simplify_visitor_0 : public postorder_visitor<expression_node,
                                                const simplify_visitor_0&,
                                                const simplify_visitor_0&>
{
  int num_of_vars;
  semantics s;
  int n, ndel;
  variable_indicator v_i;
  std::vector<unsigned int>* delnodes;
  std::vector<std::pair<unsigned int,unsigned int> >* deledges;
  unsigned int flags;
  additional_info_u m;
  int nnum;
  int l_nnum;
  unsigned int l_flags;
  semantics l_s;
  const std::vector<double> *transfer;
  model* __mod;
  int r_idx;

public:
  simplify_visitor_0() : num_of_vars(0), s(), n(0), ndel(0), v_i(0),
                        delnodes(NULL), deledges(NULL), flags(0), m(0),
                        nnum(-1), l_nnum(-1), l_flags(0), l_s(),
                        transfer(NULL), __mod(NULL), r_idx(-1) {}

public:
  simplify_visitor_0(int __nv, std::vector<unsigned int>* __dn,
                     std::vector<std::pair<unsigned int,unsigned int> >* __de,
                     model* __m) :
                        num_of_vars(__nv), s(), n(0), ndel(0), v_i(__nv),
                        delnodes(__dn), deledges(__de), flags(0), m(0),
                        nnum(-1), l_nnum(-1), l_flags(0), l_s(),
                        transfer(NULL), __mod(__m), r_idx(-1)
        { v_i.clear(); }
  // for postorder walking, the copy constructor resets the vsitor
  simplify_visitor_0(const simplify_visitor_0& __x) :
                      num_of_vars(__x.num_of_vars), s(__x.s),
                      n(__x.n), ndel(__x.ndel), v_i(__x.v_i),
                      delnodes(__x.delnodes), deledges(__x.deledges),
                      flags(__x.flags), m(__x.m), nnum(__x.nnum),
                      l_nnum(__x.l_nnum), l_flags(__x.l_flags),
                      l_s(__x.l_s), transfer(__x.transfer), __mod(__x.__mod),
                      r_idx(__x.r_idx)
        { }

  void vinit()
     { n = 0; s = semantics(); v_i.clear(); flags = 0;
       m = 0; nnum = -1; l_nnum = -1; l_flags = 0; ndel = 0; r_idx = -1;} 
  void init() { vinit(); }

    void postorder_help(const expression_node &r, unsigned int n_chld);
    bool postorder(expression_node &r);
    
    const simplify_visitor_0& value() { return *this; }
    const simplify_visitor_0& vvalue() { s.degree = -2; return *this; }

    void simple_sum_prod_update(expression_node& r, const simplify_visitor_0& __s)
    {
      if(__s.flags & SIMPLIFY_0_IS_SUM &&
         __s.flags & SIMPLIFY_0_SUM_IS_SIMPLE)
      {
        r.params = r.coeffs[0] * __s.m.nd() + r.params.nd();
        r.coeffs[0] *= (*__s.transfer)[0];
        if(__s.flags & SIMPLIFY_0_IS_GHOST)
          transfer_ghost_down(__s.nnum, r.node_num);
        else
          (*delnodes).push_back(__s.nnum);
      }
      else if(__s.flags & SIMPLIFY_0_IS_PROD &&
              __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE)
      {
        r.coeffs[0] *= __s.m.nd();
        if(__s.flags & SIMPLIFY_0_IS_GHOST)
          transfer_ghost_down(__s.nnum, r.node_num);
        else
          (*delnodes).push_back(__s.nnum);
      }
      if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
      {
        r.coeffs[n] *= __s.m.nd();
        if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
          __mod->node(__s.nnum)->params = 1.;
      }
    }

    void transfer_ghost_down(unsigned int nnum, unsigned int pnum)
    {
      // this is only called for simple sums (n+c) and simple prods (n*c)
      model::walker _w = __mod->node(nnum), _wc;
      children_iterator _c;
      std::vector<unsigned int> unum;
      for(_c = _w.child_begin(); _c != _w.child_end(); ++_c)
      {
        _wc = _w >> _c;
        unum.push_back(_wc->node_num);
      }
      // find the appropriate ghost
      model::walker _go = __mod->ghost(nnum);
      model::walker _pr = __mod->node(pnum);
      bool do_remove = _go.n_parents() == 1;
      for(unsigned int i=0; i<unum.size(); ++i)
      {
        model::walker _gn = __mod->ghost(unum[i]);
        if(i == 0)
          __mod->replace_edge_to_child(_pr, _go, _gn);
        else
          __mod->add_edge_back(_pr, _gn);
      }
      if(do_remove)
        __mod->remove_node(_go, nnum);
    }

    void vcollect(const simplify_visitor_0& __s);

    void collect(expression_node &r, const simplify_visitor_0& __s);
  };


  inline bool model::children_eq(walker __n1, walker __n2)
  {
    children_iterator _c1, _c2;
    walker _e1, _e2;

    for(_c1 = __n1.child_begin(), _c2 = __n2.child_begin();
        _c1 != __n1.child_end(); ++_c1, ++_c2)
    {
      _e1 = __n1 >> _c1;
      _e2 = __n2 >> _c2;
      if(_e1.is_sky() != _e2.is_sky())
        return false;
      if(_e1->node_num != _e2->node_num)
        return false;
    }
    return true;
  }

  inline void model::basic_simplify_1_flood(walker __s,
      std::vector<int>& _combined)
  {
    parents_iterator __p;
    walker __ew;

    for(__p = __s.parent_begin(); __p != __s.parent_end(); ++__p)
    {
      __ew = __s << __p;
      if(!__ew.is_ground())
        ++_combined[__ew->node_num];
    }
  }

  inline bool model::basic_simplify_1_is_ready(walker __s,
      const std::vector<int>& _combined, bool all_finished)
  {
    parents_iterator _p;
    walker _w;

    for(_p = __s.parent_begin(); _p != __s.parent_end(); ++_p)
    {
      _w = __s << _p;
      if(_combined[_w->node_num] == (int)_w->n_children)
      {
        if(!all_finished)
          return true;
      }
      else if(all_finished)
        return false;
    }
    return all_finished;
  }

#if 0
  class model::simplify_visitor_thin : public postorder_visitor<expression_node,
                                                                double>
  {
  private:
    double r;
    int last_nn;
    
  public:
    simplify_visitor_thin() { }
    
    void vinit() { r = INFINITY; } 

    void init() { vinit(); }

    bool postorder(expression_node &r)
    {
      if(r.operator_type < 0)
      {
        switch(r.operator_type)
        {
          case EXPRINFO_CONSTANT:
            r = r.params.nd();
            break;
          case EXPRINFO_VARIABLE:
            if(r.f_bounds.is_thin())
              r = r.f_bounds.inf();
            else
              r = INFINITY;
            break;
          case EXPRINFO_SUM:
            flags |= SIMPLIFY_0_IS_SUM;
            break;
          case EXPRINFO_PROD:
            flags |= SIMPLIFY_0_IS_PROD;
            m = r.params.nd();
            transfer = &r.coeffs;
            if(n - ndel == 1) // this is \la*<node>
              flags |= SIMPLIFY_0_PROD_IS_SIMPLE;
            break;
          case EXPRINFO_SQUARE:
          case EXPRINFO_INTPOWER:
            break;
          case EXPRINFO_MAX:
          case EXPRINFO_MIN:
          case EXPRINFO_INVERT:
          case EXPRINFO_DIV:
          case EXPRINFO_SQROOT:
          case EXPRINFO_POW:
          case EXPRINFO_ABS:
          case EXPRINFO_EXP:
          case EXPRINFO_LOG:
          case EXPRINFO_SIN:
          case EXPRINFO_COS:
          case EXPRINFO_ATAN2:
          case EXPRINFO_GAUSS:
            s.degree = -1;
            break;
          case EXPRINFO_LIN:
            break;
          case EXPRINFO_QUAD:
            break;
          default:
            // give bad messages
            break;
        }
      }
      else if(r.operator_type > 0)
      {
        s.degree = -1;
        s.property_flags.convex = t_maybe;
        s.property_flags.concave = t_maybe;
        s.property_flags.separable = t_maybe;
      }
      s.dim = v_i.sum(0, num_of_vars);
      r.sem = s;
      r.v_i = v_i;
      if(r.operator_type == EXPRINFO_CONSTANT ||
         r.operator_type == EXPRINFO_VARIABLE)
        r.n_children = 0;
      else
        r.n_children = n - ndel;
      nnum = r.node_num;
      if(s.degree < 0) s.degree = -1;
      return false;
    }
    
    const simplify_visitor_thin& value() { return *this; }
    const simplify_visitor_thin& vvalue() { s.degree = -2; return *this; }

    void simple_sum_prod_update(expression_node& r, const simplify_visitor_thin& __s)
    {
      if(__s.flags & SIMPLIFY_0_IS_SUM &&
         __s.flags & SIMPLIFY_0_SUM_IS_SIMPLE)
      {
        r.params = r.coeffs[0] * __s.m.nd();
        r.coeffs[0] *= (*__s.transfer)[0];
        (*delnodes).push_back(__s.nnum);
      }
      else if(__s.flags & SIMPLIFY_0_IS_PROD &&
              __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE)
      {
        r.coeffs[0] *= __s.m.nd();
        (*delnodes).push_back(__s.nnum);
      }
    }

    void vcollect(const simplify_visitor_thin& __s)
    {
      if((__s.flags & SIMPLIFY_0_IS_SUM &&
          __s.flags & SIMPLIFY_0_SUM_IS_SIMPLE && __s.m.nd() == 0) ||
         (__s.flags & SIMPLIFY_0_IS_PROD &&
          __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE && __s.m.nd() == 1))
      { // there is a <node>+0 or <node>*1 expression on top -> remove
        (*delnodes).push_back(__s.nnum);
      }
    }

    void collect(expression_node &r, const simplify_visitor_thin& __s)
    {
      if(r.operator_type < 0)
      {
        switch(r.operator_type)
        {
          case EXPRINFO_CONSTANT:
          case EXPRINFO_VARIABLE:
            break;
          case EXPRINFO_SUM:
            if(n == 0)
            {
              s.degree = __s.s.degree;
              s.property_flags = __s.s.property_flags;
            }
            else
            { 
              if(__s.s.degree == -1 || __s.s.degree > s.degree)
                s.degree = __s.s.degree;
            }
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { 
              double help =
                  __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                               __s.m.nd();
              r.params = r.params.nd() + r.coeffs[n-ndel] * help;
              r.coeffs.erase(r.coeffs.begin() + n);
              (*delnodes).push_back(__s.nnum);
              ndel++;
            }
            else if(__s.flags & SIMPLIFY_0_IS_PROD &&
                    __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE)
            {
              r.coeffs[n-ndel] *= __s.m.nd();
              (*delnodes).push_back(__s.nnum);
            }
            else if(__s.flags & SIMPLIFY_0_IS_SUM)
            {
              double help = r.coeffs[n-ndel];

              r.params = r.params.nd() + __s.m.nd();
              r.coeffs.erase(r.coeffs.begin()+(n-ndel));
              ndel++;
              for(std::vector<double>::iterator b = __s.transfer->begin();
                  b != __s.transfer->end(); ++b)
                r.coeffs.push_back(*b * help);
              n += __s.transfer->size();
              (*delnodes).push_back(__s.nnum);
            }
            break;
          case EXPRINFO_PROD:
            if(n == 0)
            {
              s.degree = __s.s.degree;
              s.property_flags = __s.s.property_flags;
            }
            else
            { 
              if(__s.s.degree == -1 || s.degree == -1)
                s.degree = -1;
              else
                s.degree += __s.s.degree;
            }
            if(r.coeffs[n] != 1)
            {
              r.params = r.params.nd() * r.coeffs[n];
              r.coeffs[n] = 1;
            }
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { 
              double help =
                  __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                               __s.m.nd();
              r.params = r.params.nd() * help;
              r.coeffs.pop_back();
              (*delnodes).push_back(__s.nnum);
              ndel++;
            }
            else if(__s.flags & SIMPLIFY_0_IS_PROD)
            {
              r.params = r.params.nd() * __s.m.nd();
              r.coeffs.insert(r.coeffs.end(), __s.transfer->size()-1, 1.0);
              n += __s.transfer->size()-1;
              (*delnodes).push_back(__s.nnum);
            }
            break;
          case EXPRINFO_MAX:
          case EXPRINFO_MIN:
            break;
          case EXPRINFO_POW:
            if(n == 0)
            {
              if(__s.flags & SIMPLIFY_0_IS_CONST)
              { // a^x = exp(x * ln a)
                double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                                 __s.m.nd();
                r.operator_type = EXPRINFO_EXP;
                r.coeffs[0] = log(r.coeffs[0] * help) * r.coeffs[1];
                r.coeffs.pop_back();
                (*delnodes).push_back(__s.nnum);
                s.degree = -1;
                ndel++;
              }
              else
                s.degree = __s.s.degree;
              simple_sum_prod_update(r, __s);
            }
            if(n != 0 && r.params.nd() == 0 &&
               __s.flags & SIMPLIFY_0_IS_CONST &&
               __s.flags & SIMPLIFY_0_CONST_IS_INTEGER)
            {
              r.params = __s.m.nn();
              r.operator_type = EXPRINFO_INTPOWER;
              (*delnodes).push_back(__s.nnum);
              s.degree = s.degree == -1 ? -1 : s.degree*__s.m.nn();
              ndel++;
            }
            break;
          case EXPRINFO_SQUARE:
            s.degree = __s.s.degree == -1 ? -1 : __s.s.degree*2;
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help += r.params.nd();
              r.params = help*help;
              (*delnodes).push_back(__s.nnum);
            }
            else
              simple_sum_prod_update(r, __s);
            break;
          case EXPRINFO_SQROOT:
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help += r.params.nd();
              r.params = sqrt(help);
              (*delnodes).push_back(__s.nnum);
            }
            else
              simple_sum_prod_update(r, __s);
            break;
          case EXPRINFO_ABS:
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help += r.params.nd();
              r.params = fabs(help);
              (*delnodes).push_back(__s.nnum);
            }
            else
              simple_sum_prod_update(r, __s);
            break;
          case EXPRINFO_EXP:
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help += r.params.nd();
              r.params = exp(help);
              (*delnodes).push_back(__s.nnum);
            }
            else
              simple_sum_prod_update(r, __s);
            break;
          case EXPRINFO_LOG:
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help += r.params.nd();
              r.params = log(help);
              (*delnodes).push_back(__s.nnum);
            }
            else
              simple_sum_prod_update(r, __s);
            break;
          case EXPRINFO_SIN:
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help += r.params.nd();
              r.params = sin(help);
              (*delnodes).push_back(__s.nnum);
            }
            else
              simple_sum_prod_update(r, __s);
            break;
          case EXPRINFO_COS:
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help += r.params.nd();
              r.params = cos(help);
              (*delnodes).push_back(__s.nnum);
            }
            else
              simple_sum_prod_update(r, __s);
            break;
          case EXPRINFO_ATAN2:
            // const should be computed here, too.
            break;
          case EXPRINFO_GAUSS:
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              help *= r.coeffs[0];
              help -= r.params.d()[0];
              help /= r.params.d()[1];
              r.params = exp(-help*help);
              (*delnodes).push_back(__s.nnum);
            }
            else if(__s.flags & SIMPLIFY_0_IS_SUM &&
               __s.flags & SIMPLIFY_0_SUM_IS_SIMPLE)
            {
              std::vector<double> m = r.params.d();

              m[0] -= r.coeffs[0]*__s.m.nd();
              r.params = m;
              r.coeffs[0] *= (*__s.transfer)[0];
              (*delnodes).push_back(__s.nnum);
            }
            else if(__s.flags & SIMPLIFY_0_IS_PROD &&
                    __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE)
            {
              r.coeffs[0] *= __s.m.nd();
              (*delnodes).push_back(__s.nnum);
            }
            break;
          case EXPRINFO_INVERT:
            if(r.coeffs[0] != 1)
            {
              r.params = r.params.nd() / r.coeffs[0];
              r.coeffs[0] = 1.;
            }
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              r.params = r.params.nd() / help;
              (*delnodes).push_back(__s.nnum);
            }
            break;
          case EXPRINFO_DIV:
            if(n == 0 && __s.flags & SIMPLIFY_0_IS_CONST)
            {
              double help =
                  __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                               __s.m.nd();
              r.operator_type = EXPRINFO_INVERT;
              r.params = help*r.coeffs[0];
              r.coeffs.erase(r.coeffs.begin());
              (*delnodes).push_back(__s.nnum);
              s.degree = __s.s.degree;
              ndel++;
            }
            else if(n == 0)
            {
              r.params = r.coeffs[0];
              r.coeffs[0] = 1;
            }
            if(n != 0)
            {
              if(__s.flags & SIMPLIFY_0_IS_CONST)
              {
                double help =
                     __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                               __s.m.nd();
                r.operator_type = EXPRINFO_PROD;
                r.params = 1./help;
                r.coeffs.pop_back();
                (*delnodes).push_back(__s.nnum);
                ndel++;
              }
              else
              {
                r.params = r.params.nd() / r.coeffs[1];
                r.coeffs[1] = 1;
                s.degree = -1;
              }
            }
            break;
          case EXPRINFO_INTPOWER:
            if(r.params.nn() > 0)
              s.degree = __s.s.degree == -1 ? -1 : __s.s.degree*r.params.nn();
            else
              s.degree = -1;
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // compute constant expressions
              double help =
                    __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                              __s.m.nd();
              r.operator_type = EXPRINFO_CONSTANT;
              r.params = pow(help, r.params.nn());
              (*delnodes).push_back(__s.nnum);
            }
            break;
          case EXPRINFO_LIN:
            break;
          case EXPRINFO_QUAD:
            break;
          default:
            // give bad messages
            break;
        }
      }
      else if(r.operator_type > 0)
      {
        ;
      }
      n++;
      v_i.set_all(__s.v_i);
    }
  };
#endif

  inline bool model::simplify_thin()
  {
    return false;
  }

  inline void model::arrange_constraints()
  {
    sort(constraints.begin(), constraints.end(), sort_constraints());
  }

  inline int model::next_num() { return gid->get_node_id(); }
  inline int model::next_variable_num() { return gid->get_var_id(); }
  inline int model::next_constraint_num() { return gid->get_const_id(); }

  inline void model::renumber_variables()
  {
    gid->compress_numbers(true);
    detect_0chain();
  }

  inline void model::compress_numbers()
  {
    gid->compress_numbers();
    detect_0chain();
  }

  struct model::__docompare_nodes
  {
    bool operator() (const model::walker& _a, unsigned int _b) const
      {return _a->node_num < _b;}
    bool operator() (unsigned int _a, const model::walker& _b) const
      {return _a < _b->node_num;}
    bool operator() (const model::walker& _a, const model::walker& _b) const
      {return _a->node_num < _b->node_num;}
  };

  struct model::__docompare_variables
  {
    bool operator() (const model::walker& _a, unsigned int _b) const
      {return _a->params.nn() < (int)_b;}
    bool operator() (unsigned int _a, const model::walker& _b) const
      {return (int)_a < _b->params.nn();}
    bool operator() (const model::walker& _a, const model::walker& _b) const
      {return _a->params.nn() < _b->params.nn();}
  };

  inline model::walker model::store_node(const walker& _w)
  {
    if(!gid->have_glob_ref(_w->node_num))
      gid->number_of_nodes(_w->node_num+1);

    ref_iterator __x(lower_bound(node_ref.begin(),
                  node_ref.end(), _w->node_num, __docompare_nodes()));

    if(__x != node_ref.end() && (*__x)->node_num == _w->node_num)
      *__x = _w;
    else
      node_ref.insert(__x, _w);
    gid->mk_globref(_w->node_num, _w);
    return _w;
  }

  inline model::walker model::store_variable(const walker& _w)
  {
    if(_w->node_num >= number_of_nodes())
      gid->number_of_nodes(_w->node_num+1);

    ref_iterator __x(lower_bound(var_ref.begin(),
                  var_ref.end(), _w->params.nn(), __docompare_variables()));

    if(__x != var_ref.end() && (*__x)->params.nn() == _w->params.nn())
      *__x = _w;
    else
      var_ref.insert(__x, _w);
    gid->mk_gvarref(_w->params.nn(), _w);
    return _w;
  }

  inline model::walker model::store_constraint(const walker& _w)
  {
    unsigned int dummy;
    if(!gid->get_const_num(_w->node_num, dummy))
    {
      // it is not already stored
      constraints.push_back(_w);
      gid->mk_gconstref(next_constraint_num(), _w);
    }
    return _w;
  }

  inline model::walker model::store_ghost(const walker& _w)
  {
    if(!gid->have_glob_ref(_w->node_num))
      gid->number_of_nodes(_w->node_num+1);

    ref_iterator __x(lower_bound(ghost_ref.begin(),
                  ghost_ref.end(), _w->node_num, __docompare_nodes()));

    if(__x != ghost_ref.end() && (*__x)->node_num == _w->node_num)
      *__x = _w;
    else
      ghost_ref.insert(__x, _w);
    return _w;
  }

  inline void model::free_node_num(unsigned int _nnum)
  {
#if MODEL_INLINE_DEBUG
    std::cout << "Removing node " << _nnum << std::endl;
#endif
    ref_iterator __x(lower_bound(node_ref.begin(),
                  node_ref.end(), _nnum, __docompare_nodes()));
    if(__x != node_ref.end() && (*__x)->node_num == _nnum)
      node_ref.erase(__x);
    gid->remove_node_id(_nnum);
  }

  inline void model::remove_node(const walker& _w, unsigned int _nnum)
  {
    walker __h = _w;
    unsigned int cn;

    if(_w->operator_type == EXPRINFO_GHOST)
    { // only remove ghosts from the reference
      ref_iterator __x(lower_bound(ghost_ref.begin(),
                    ghost_ref.end(), _nnum, __docompare_nodes()));
      if(__x != ghost_ref.end() && (*__x)->node_num == _nnum)
        ghost_ref.erase(__x);
      erase(__h);
      return;
    }
    if(is_empty(_w))
    {
      std::cerr << "Error: Trying to remove node, which is already removed.";
      std::terminate();
      return;
    }
    if(_w == objective)
    {
      if(_w->n_children != 1)
      // cannot erase the objective node unless it has only one child
      // never mind about not removing it, it has to be a root node
      // anyway.
        return;
      walker _ew = _w;
      objective = _ew >> _ew.child_begin();
    }
    if(get_const_num(_nnum, cn))
    {
      unsigned int ccn;
      if(_w->n_children != 0)
      {
        walker _ew = _w;
        _ew >>= _ew.child_begin();
        if(get_const_num(_ew->node_num, ccn)) // the child is a constraint as well
        {
          gid->unused_constr(const_name(cn));
          gid->remove_const_id(cn);
        }
        else
          gid->make_const_back_ref(_ew->node_num, cn);
      }
      else
      {
        gid->unused_constr(const_name(cn));
        gid->remove_const_id(cn);
      }
      // remove the constraint from the constraints array
      for(std::vector<walker>::iterator __x = constraints.begin();
          __x != constraints.end(); ++__x)
        if(*__x == _w)
        {
          constraints.erase(__x);
          break;
        }
    }
    ref_iterator __x(lower_bound(node_ref.begin(),
                  node_ref.end(), _nnum, __docompare_nodes()));
    if(__x != node_ref.end() && (*__x)->node_num == _nnum)
      node_ref.erase(__x);
    free_node_num(_nnum);
    erase(__h);
  }

  inline void model::remove_node(const walker& _w)
  {
    remove_node(_w, _w->node_num);
  }

  inline void model::remove_node(unsigned int __node_num)
  {
    remove_node(gid->node(__node_num), __node_num);
  }

  inline void model::new_variables(int _new_num_of_vars)
  {
    gid->number_of_variables(_new_num_of_vars);
  }

  inline model::walker model::constant(double _constant)
  {
    expression_node _en(EXPRINFO_CONSTANT, next_num());

    _en.v_i.reserve(gid->number_of_variables());
    _en.v_i.clear();
    _en.params = _constant;

    return store_node(between_back(ground(), sky(), _en));
  }

  inline model::walker model::constant(const std::vector<double>& _constant)
  {
    expression_node _en(EXPRINFO_CONSTANT, next_num());

    _en.v_i.reserve(gid->number_of_variables());
    _en.v_i.clear();
    _en.params = _constant;

    return store_node(between_back(ground(), sky(), _en));
  }

  inline model::walker model::ghost(unsigned int _nnum)
  {
    ref_iterator __x;

    __x = lower_bound(ghost_ref.begin(), ghost_ref.end(), _nnum,
                        __docompare_nodes());
    if(__x != ghost_ref.end() && (*__x)->node_num == _nnum)
      return *__x;
    else // make a new one
    {
      walker __r;
      expression_node _en(EXPRINFO_GHOST, _nnum);
      _en.params = false;
      _en.f_bounds = interval(-INFINITY,INFINITY);
      __r = between_back(ground(), sky(), _en);
      ghost_ref.insert(__x, __r);
      return __r;
    }
  }

  inline model::walker model::binary(const walker& _op1, const walker& _op2,
                              int expr_type, double _coeff1, double _coeff2)
  {
    expression_node _en(expr_type, next_num());
    std::vector<walker> _w(2);

    _en.v_i = (*_op1).v_i;
    _en.v_i.set_all((*_op2).v_i);
    _w[0] = _op1;
    _w[1] = _op2;
    _en.coeffs.push_back(_coeff1);
    _en.coeffs.push_back(_coeff2);
    if(expr_type == EXPRINFO_PROD || expr_type == EXPRINFO_DIV||
       expr_type == EXPRINFO_INVERT)
      _en.params = 1.;
    else
      _en.params = 0.;
    return store_node(split_back(ground(), _w, _en));
  }

  inline model::walker model::binary(const walker& _op1, const walker& _op2,
                              int expr_type, additional_info_u _params,
                              double _coeff1, double _coeff2)
  {
    expression_node _en(expr_type, next_num());
    std::vector<walker> _w(2);

    _en.v_i = (*_op1).v_i;
    _en.v_i.set_all((*_op2).v_i);
    _en.params = _params;
    _w[0] = _op1;
    _w[1] = _op2;
    _en.coeffs.push_back(_coeff1);
    _en.coeffs.push_back(_coeff2);
    return store_node(split_back(ground(), _w, _en));
  }

  inline model::walker model::unary(const walker& _op1, int expr_type,
                             double _coeff)
  {
    expression_node _en(expr_type, next_num());

    _en.v_i = (*_op1).v_i;
    _en.coeffs.push_back(_coeff);
    if(expr_type == EXPRINFO_INVERT)
      _en.params = 1.;
    else
      _en.params = 0.;
    return store_node(split_back(ground(), _op1, _en));
  }

  inline model::walker model::unary(const walker& _op1, int expr_type,
                             additional_info_u _params, double _coeff)
  {
    expression_node _en(expr_type, next_num());

    _en.v_i = (*_op1).v_i;
    _en.params = _params;
    _en.coeffs.push_back(_coeff);
    return store_node(split_back(ground(), _op1, _en));
  }

  inline model::walker model::nary(const std::vector<walker>& _op, int expr_type,
                                   const std::vector<double>& _coeffs)
  {
    expression_node _en(expr_type, next_num());

    _en.v_i = (*_op[0]).v_i;
    for(const_ref_iterator b = _op.begin(); b != _op.end(); ++b)
      _en.v_i.set_all((**b).v_i);
    if(_coeffs.size() == 0)
      _en.coeffs.insert(_en.coeffs.end(),_op.size(),1.0);
    else
      _en.coeffs = _coeffs;
    if(expr_type == EXPRINFO_PROD)
      _en.params = 1.;
    else
      _en.params = 0.;
    return store_node(split_back(ground(), _op, _en));
  }

  inline model::walker model::nary(const std::vector<walker>& _op, int expr_type,
                            additional_info_u _params,
                            const std::vector<double>& _coeffs)
  {
    expression_node _en(expr_type, next_num());

    _en.v_i = (*_op[0]).v_i;
    for(const_ref_iterator b = _op.begin(); b != _op.end(); ++b)
      _en.v_i.set_all((**b).v_i);
    _en.params = _params;
    if(_coeffs.size() == 0)
      _en.coeffs.insert(_en.coeffs.end(),_op.size(),1.0);
    else
      _en.coeffs = _coeffs;
    return store_node(split_back(ground(), _op, _en));
  }

  inline void model::clr_sky_ground_link()
  {
    if(ground().n_children() != 1 || !empty())
      remove_edge(ground(), sky());
  }

  inline void model::set_counters()
  {
    ref_iterator __e = node_ref.end();
    walker _x;
    for(ref_iterator __x = node_ref.begin(); __x != __e; ++__x)
    {
      _x = *__x;
      if(_x.n_children() == 1)
      {
        walker _c = _x >> _x.child_begin();
        if(_c == sky())
          _x->n_children = 0;
        else
          _x->n_children = 1;
      }
      else
        _x->n_children = _x.n_children();
      if(_x.n_parents() == 1)
      {
        walker _c = _x >> _x.parent_begin();
        if(_c == ground())
          _x->n_parents = 0;
        else
          _x->n_parents = 1;
      }
      else
        _x->n_parents = _x.n_parents();
    }
  }

  struct model::detect_0chain_visitor_st
  {
    std::pair<unsigned int,unsigned int> r;    // return value (node num, depth)
    unsigned int *f;                      // counter
    std::vector<unsigned int>* t;              // type
    std::vector<unsigned int>* d;              // depth
    std::vector<unsigned int>* nn;             // base
    bool new_number;                      // a new chain??
  };

  class model::detect_0chain_visitor : public
    cached_forward_evaluator_base<model::detect_0chain_visitor_st,
                          expression_node, std::pair<unsigned int,unsigned int>,
                          model::walker>
  {
  private:
    typedef cached_forward_evaluator_base<model::detect_0chain_visitor_st,
                  expression_node, std::pair<unsigned int,unsigned int>,
                  model::walker> _Base;
  public:
    detect_0chain_visitor(std::vector<unsigned int>& __t, std::vector<unsigned int>& __d,
                             std::vector<unsigned int>& __nn, unsigned int &d)
          {
            eval_data.r = std::make_pair(0,0);
            eval_data.t = &__t;
            eval_data.d = &__d;
            eval_data.f = &d;
            eval_data.nn = &__nn;
            eval_data.new_number = false;
          }

    detect_0chain_visitor(const detect_0chain_visitor& __x) : _Base(__x)
          {} // constructor os assign the value

    ~detect_0chain_visitor() {} // destructor

    void initialize() { return; } // does nothing

    bool is_cached(const expression_node& __data)
      { return (*eval_data.t)[__data.node_num] != 0; }

    void retrieve_from_cache(const expression_node& __data)
      { eval_data.r = std::make_pair((*eval_data.t)[__data.node_num],
                                (*eval_data.d)[__data.node_num]); }

    int initialize(const expression_node& __data)
      { 
        if(__data.operator_type == EXPRINFO_VARIABLE)
        {
          unsigned int __h;
          eval_data.new_number = true;
          __h = (*eval_data.f)++;
          eval_data.r = std::make_pair(__h,0);
          (*eval_data.nn).push_back(__data.node_num);
        }
        else
        {
          eval_data.new_number = false;
          eval_data.r = std::make_pair(0,0);
        }
        return 1;
      }

    void calculate(const expression_node& __data)
      {
        (*eval_data.t)[__data.node_num] = eval_data.r.first;
        (*eval_data.d)[__data.node_num] = eval_data.r.second;
      }

    int update(const std::pair<unsigned int,unsigned int>& __rval) { return 0; }

    int update(const expression_node& __data,
                const std::pair<unsigned int,unsigned int>& __rval)
      {
        if(__data.operator_type == EXPRINFO_CONSTANT ||
           __data.operator_type == EXPRINFO_VARIABLE || __rval.first == 0)
          return 0;
        if(eval_data.r.first == 0)
          eval_data.r = std::make_pair(__rval.first,__rval.second+1);
        else if(eval_data.r.first == __rval.first)
          eval_data.r.second++;  // so it definitely is >1 later.
        else if(!eval_data.new_number)
        {
          int __h;
          eval_data.new_number = true;
          __h = (*eval_data.f)++;
          eval_data.r = std::make_pair(__h,0);
          (*eval_data.nn).push_back(__data.node_num);
        }
        return 0;
      }

    std::pair<unsigned int,unsigned int> calculate_value(bool eval_all)
      { return eval_data.r; }
  };

  struct model::lincoeff_visitor_ret
  {
    bool was_lin;
    bool was_const;
    int node_num;
    double const_trans;
  };

  struct model::lincoeff_visitor_st
  {
    const model* mod;
    const std::vector<interval>* ranges;
    sparse_vector<double>* coeffs;
    //mtl::compressed1D<double>* coeffs;
    double* constant;
    double old_trans, trans, mm_help;
    lincoeff_visitor_ret r;
    int old_nlin, in_nlin;
    unsigned int n;
  };

  class model::lincoeff_visitor : public
    cached_backward_evaluator_base<lincoeff_visitor_st, expression_node,
                            lincoeff_visitor_ret, model::walker>
  {
  private:
    typedef cached_backward_evaluator_base<lincoeff_visitor_st, expression_node,
                              lincoeff_visitor_ret, model::walker> _Base;

  public:
    lincoeff_visitor(const std::vector<interval>* _rg,
                     sparse_vector<double>* _cf, double* _cs,
                     const model* _m)
    {
      eval_data.trans = 1;
      eval_data.old_trans = 1;
      eval_data.n = 0;
      eval_data.r.was_lin = true;
      eval_data.r.was_const = false;
      eval_data.r.const_trans = 0;
      eval_data.r.node_num = -1;
      eval_data.old_nlin = -1;
      eval_data.in_nlin = -1;
      eval_data.ranges = _rg;
      eval_data.coeffs = _cf;
      eval_data.constant = _cs;
      eval_data.mod = _m;
    }

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

    ~lincoeff_visitor() {}

    bool is_cached(const expression_node& __data)
    {
      return (__data.operator_type == EXPRINFO_LIN ||
          __data.operator_type == EXPRINFO_QUAD);
    }

    void initialize()
    {
      eval_data.n = 0;
      eval_data.r.was_lin = false;
      eval_data.r.was_const = false;
      eval_data.r.node_num = -1;
      eval_data.r.const_trans = 0;
      eval_data.old_nlin = eval_data.in_nlin;
      eval_data.old_trans = eval_data.trans;
    }

    void initialize(const expression_node& __data)
    {
      eval_data.n = 0;
      eval_data.r.was_lin = false;
      eval_data.r.was_const = false;
      eval_data.r.node_num = __data.node_num;
      eval_data.r.const_trans = 0;
      eval_data.old_nlin = eval_data.in_nlin;
      eval_data.old_trans = eval_data.trans;
    }

    void retrieve_from_cache(const expression_node& __data)
    {
      if(__data.operator_type == EXPRINFO_LIN)
      {
        eval_data.r.was_lin = true;
        if(eval_data.trans == 1)
        {
          linalg_add(*eval_data.coeffs, eval_data.mod->lin[__data.params.nn()],
              *eval_data.coeffs);
        }
        else
        {
          linalg_add(*eval_data.coeffs,
              linalg_scale(eval_data.mod->lin[__data.params.nn()],
                           eval_data.trans),
              *eval_data.coeffs);
        }
      }
      else if(__data.operator_type == EXPRINFO_QUAD)
        eval_data.r.was_lin = eval_data.r.was_const = false;
    }

    int calculate(const expression_node& __data)
    {
      if(__data.operator_type < 0)
      {
        switch(__data.operator_type)
        {
          case EXPRINFO_CONSTANT:
            eval_data.r.was_lin = true;
            eval_data.r.was_const = true;
            eval_data.r.const_trans = __data.params.nd();
            break;
          case EXPRINFO_VARIABLE:
            if(eval_data.ranges &&
               (*eval_data.ranges)[__data.node_num].is_thin())
            {
              eval_data.r.was_const = true;
              eval_data.r.const_trans =
                                  (*eval_data.ranges)[__data.node_num].sup();
            }
            else if(__data.f_bounds.is_thin())
            {
              eval_data.r.was_const = true;
              eval_data.r.const_trans = __data.f_bounds.sup();
            }
            else if(eval_data.old_nlin == -1)
              (*eval_data.coeffs)[__data.params.nn()] += eval_data.old_trans;
            eval_data.r.was_lin = true;
            break;
          case EXPRINFO_SUM:
            if(eval_data.old_nlin == -1)
              *eval_data.constant += eval_data.old_trans * __data.params.nd();
            break;
          case EXPRINFO_PROD:
          case EXPRINFO_DIV:
            eval_data.old_trans *= __data.params.nd();
            eval_data.in_nlin = 1;
            break;
          case EXPRINFO_INVERT:
            eval_data.old_trans *= __data.params.nd();
            eval_data.in_nlin = 2;
            break;
          case EXPRINFO_MAX:
          case EXPRINFO_MIN:
            eval_data.in_nlin = -2;
            eval_data.r.was_lin = eval_data.r.was_const = true;
            eval_data.r.const_trans = __data.params.nd();
            break;
          case EXPRINFO_ABS:
            if(eval_data.ranges)
            {
              if((*eval_data.ranges)[__data.node_num].sup() <= 0)
                eval_data.old_trans *= -1;
              else if((*eval_data.ranges)[__data.node_num].inf() < 0)
                eval_data.in_nlin = 1;
            }
            else
            {
              if(__data.f_bounds.sup() <= 0)
              {
                eval_data.old_trans *= -1;
              }
              else if(__data.f_bounds.inf() < 0)
                eval_data.in_nlin = 1;
            }
            break;
          case EXPRINFO_ATAN2:
            eval_data.in_nlin = 2;
            break;
          case EXPRINFO_POW:
            eval_data.in_nlin = 2;
            break;
          case EXPRINFO_INTPOWER:
            if(__data.params.nn() != 1)
              eval_data.in_nlin = 2;
            break;
          case EXPRINFO_SQUARE:
          case EXPRINFO_SQROOT:
          case EXPRINFO_EXP:
          case EXPRINFO_LOG:
          case EXPRINFO_SIN:
          case EXPRINFO_COS:
          case EXPRINFO_GAUSS:
            eval_data.in_nlin = 2;
            break;
          case EXPRINFO_LIN:
          case EXPRINFO_QUAD:
            std::cerr << "Must not be here in lincoeff_visitor" << std::endl;
            throw "Strange error";
            break;
          default:
            std::cerr << "Unknown expression type in get_linear_coeffs" <<
              std::endl;
            throw "Programming error";
            break;
        }
      }
      else if(__data.operator_type > 0)
      {
        eval_data.r.was_lin = false;
        eval_data.r.was_const = false;
        eval_data.in_nlin = 3;
      }
      if(!__data.coeffs.empty())
        eval_data.trans = eval_data.old_trans * __data.coeffs[0];
      return 1;
    }

    int update(const expression_node& __data, const lincoeff_visitor_ret& __rval)
    {
      if(__data.operator_type < 0)
      {
        switch(__data.operator_type)
        {
          case EXPRINFO_CONSTANT:
          case EXPRINFO_VARIABLE:
            break;
          case EXPRINFO_SUM:
            if(eval_data.n == 0)
            {
              eval_data.r.was_lin = __rval.was_lin;
              eval_data.r.was_const = __rval.was_const;
            }
            else if(__rval.was_const && eval_data.r.was_const)
            {
              eval_data.r.was_lin = true;
              eval_data.r.const_trans = 
                (eval_data.n == 0 ? __rval.const_trans
                           : __rval.const_trans + eval_data.r.const_trans);
            }
            else if(__rval.was_const && eval_data.r.was_lin)
              eval_data.r.const_trans += __rval.const_trans;
            else if(eval_data.r.was_const && __rval.was_lin)
            {
              eval_data.r.was_const = false;
              eval_data.r.was_lin = true;
            }
            else if(!__rval.was_lin)
              eval_data.r.was_lin = false;
            if(eval_data.n == __data.n_children - 1)
            {
              if(eval_data.r.was_const)
                eval_data.r.const_trans += __data.params.nd();
              else if(eval_data.r.was_lin && eval_data.old_nlin == -1)
                *eval_data.constant += eval_data.old_trans*
                  (eval_data.r.const_trans + __data.params.nd());
            }
            break;
          case EXPRINFO_PROD:
            if(eval_data.n == 0)
              eval_data.r = __rval;
            else if(__rval.was_const && eval_data.r.was_const)
            {
              eval_data.r.was_lin = true;
              eval_data.r.const_trans = 
                (eval_data.n == 0 ? __rval.const_trans
                           : __rval.const_trans * eval_data.r.const_trans);
            }
            else if(__rval.was_const && eval_data.r.was_lin)
              eval_data.r.const_trans *= __rval.const_trans;
            else if(eval_data.r.was_const && __rval.was_lin)
            {
              eval_data.r.was_const = false;
              eval_data.r.was_lin = true;
              eval_data.r.node_num = __rval.node_num;
            }
            else
              eval_data.r.was_const = eval_data.r.was_lin = false;
            if(eval_data.n == __data.n_children - 1)
            {
              if(eval_data.r.was_lin && !eval_data.r.was_const)
              {
                lincoeff_visitor lv(eval_data.ranges, eval_data.coeffs,
                                    eval_data.constant, eval_data.mod);
                lv.eval_data.trans =
                                  eval_data.old_trans * eval_data.r.const_trans;
                eval_data.r = evaluate(lv,
                                  eval_data.mod->node(eval_data.r.node_num));
              }
              else if(eval_data.r.was_const)
                eval_data.r.const_trans *= __data.params.nd();
              eval_data.r.node_num = __data.node_num;
            }
            break;
          case EXPRINFO_INVERT:
            if(__rval.was_const)
            {
              eval_data.r.was_const = true;
              eval_data.r.was_lin = true;
              eval_data.r.const_trans =
                          __data.params.nd()/__rval.const_trans;
            }
            else
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_DIV:
            if(__rval.was_const)
            {
              if(eval_data.n == 0)
              {
                eval_data.r.const_trans = __rval.const_trans;
                eval_data.r.was_const = true;
              }
              else if(eval_data.n == 1)
              {
                if(eval_data.r.was_const && __rval.const_trans != 0)
                {
                  eval_data.r.const_trans =
                          eval_data.r.const_trans/__rval.const_trans;
                  eval_data.r.was_lin = true;
                }
                else
                  eval_data.r.was_lin = eval_data.r.was_const = false;
              }
            }
            else
              eval_data.r.was_const = eval_data.r.was_lin = false;
            break;
          case EXPRINFO_MAX:
            if(__rval.was_const && eval_data.r.was_const)
            {
              if(__rval.const_trans > eval_data.r.const_trans)
                eval_data.r.const_trans = __rval.const_trans;
            }
            else if(__rval.was_const && !eval_data.r.was_const)
            {
              if(__rval.const_trans > eval_data.mm_help)
              {
                eval_data.r.was_const = true;
                eval_data.r.was_lin = true;
                eval_data.r.const_trans = __rval.const_trans;
              }
            }
            else if(eval_data.r.was_const && !__rval.was_const)
            {
              interval _rg =
                 (*eval_data.ranges)[__rval.node_num]*__data.coeffs[eval_data.n];
              if(eval_data.r.const_trans <= _rg.sup())
              {
                eval_data.r.was_const = false;
                eval_data.r.was_lin = __rval.was_lin;
                eval_data.r.node_num = __rval.node_num;
                eval_data.r.const_trans = _rg.inf();
                eval_data.mm_help = _rg.sup();
                eval_data.in_nlin = eval_data.n;
              }
            }
            else // both were not const
            {
              interval _rg =
                 (*eval_data.ranges)[__rval.node_num]*__data.coeffs[eval_data.n];
              if(eval_data.r.const_trans <= _rg.sup()) // otherwise no change
              {
                if(eval_data.mm_help < _rg.inf()) // the new one rules
                {
                  eval_data.r.was_lin = __rval.was_lin;
                  eval_data.r.node_num = __rval.node_num;
                  eval_data.r.const_trans = _rg.inf();
                  eval_data.mm_help = _rg.sup();
                  eval_data.in_nlin = eval_data.n;
                }
                else // they mix
                {
                  eval_data.r.was_lin = false;
                  if(_rg.inf() < eval_data.r.const_trans)
                    eval_data.r.const_trans = _rg.inf();
                  if(_rg.sup() > eval_data.mm_help)
                    eval_data.mm_help = _rg.sup();
                }
              }
            }
            if(eval_data.n == __data.n_children - 1)
            {
              if(eval_data.r.was_lin && !eval_data.r.was_const &&
                 eval_data.old_nlin == -1)
              {
                lincoeff_visitor lv(eval_data.ranges, eval_data.coeffs,
                                    eval_data.constant, eval_data.mod);
                lv.eval_data.trans =
                          eval_data.old_trans * __data.coeffs[eval_data.in_nlin];
                eval_data.r = evaluate(lv,
                                  eval_data.mod->node(eval_data.r.node_num));
              }
              eval_data.r.node_num = __data.node_num;
            }
            break;
          case EXPRINFO_MIN:
            if(__rval.was_const && eval_data.r.was_const)
            {
              if(__rval.const_trans < eval_data.r.const_trans)
                eval_data.r.const_trans = __rval.const_trans;
            }
            else if(__rval.was_const && !eval_data.r.was_const)
            {
              if(__rval.const_trans < eval_data.mm_help)
              {
                eval_data.r.was_const = true;
                eval_data.r.was_lin = true;
                eval_data.r.const_trans = __rval.const_trans;
              }
            }
            else if(eval_data.r.was_const && !__rval.was_const)
            {
              interval _rg =
                 (*eval_data.ranges)[__rval.node_num]*__data.coeffs[eval_data.n];
              if(eval_data.r.const_trans >= _rg.inf())
              {
                eval_data.r.was_const = false;
                eval_data.r.was_lin = __rval.was_lin;
                eval_data.r.node_num = __rval.node_num;
                eval_data.r.const_trans = _rg.sup();
                eval_data.mm_help = _rg.inf();
                eval_data.in_nlin = eval_data.n;
              }
            }
            else // both were not const
            {
              interval _rg =
                 (*eval_data.ranges)[__rval.node_num]*__data.coeffs[eval_data.n];
              if(eval_data.r.const_trans >= _rg.inf()) // otherwise no change
              {
                if(eval_data.mm_help > _rg.sup()) // the new one rules
                {
                  eval_data.r.was_lin = __rval.was_lin;
                  eval_data.r.node_num = __rval.node_num;
                  eval_data.r.const_trans = _rg.sup();
                  eval_data.mm_help = _rg.inf();
                  eval_data.in_nlin = eval_data.n;
                }
                else // they mix
                {
                  eval_data.r.was_lin = false;
                  if(_rg.sup() > eval_data.r.const_trans)
                    eval_data.r.const_trans = _rg.sup();
                  if(_rg.inf() < eval_data.mm_help)
                    eval_data.mm_help = _rg.inf();
                }
              }
            }
            if(eval_data.n == __data.n_children - 1)
            {
              if(eval_data.r.was_lin && !eval_data.r.was_const &&
                 eval_data.old_nlin == -1)
              {
                lincoeff_visitor lv(eval_data.ranges, eval_data.coeffs,
                                    eval_data.constant, eval_data.mod);
                lv.eval_data.trans =
                          eval_data.old_trans * __data.coeffs[eval_data.in_nlin];
                eval_data.r = evaluate(lv,
                          eval_data.mod->node(eval_data.r.node_num));
              }
              eval_data.r.node_num = __data.node_num;
            }
            break;
          case EXPRINFO_ABS:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = fabs(__data.coeffs[0]*__rval.const_trans+
                  __data.params.nd());
            }
            else if(eval_data.ranges && eval_data.old_nlin == -1)
            {
              if((*eval_data.ranges)[__data.node_num].sup() <= 0)
                *eval_data.constant -= eval_data.old_trans*__data.params.nd();
              else if((*eval_data.ranges)[__data.node_num].inf() >= 0)
                *eval_data.constant += eval_data.old_trans*__data.params.nd();
            }
            else if(eval_data.old_nlin == -1)
            {
              if(__data.f_bounds.sup() <= 0)
                *eval_data.constant -= eval_data.old_trans*__data.params.nd();
              else if(__data.f_bounds.inf() >= 0)
                *eval_data.constant += eval_data.old_trans*__data.params.nd();
            }
            break;
          case EXPRINFO_SQUARE:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = __data.coeffs[0]*__rval.const_trans+
                  __data.params.nd();
              eval_data.r.const_trans *= eval_data.r.const_trans;
            }
            else
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_SQROOT:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = sqrt(__data.coeffs[0]*__rval.const_trans+
                  __data.params.nd());
            }
            else
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_INTPOWER:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              switch(__data.params.nn())
              {
                case 0:
                  eval_data.r.const_trans = 1;
                  break;
                case 1:
                  eval_data.r.const_trans = __data.coeffs[0]*__rval.const_trans;
                  break;
                case 2:
                  eval_data.r.const_trans = __data.coeffs[0]*__rval.const_trans;
                  eval_data.r.const_trans *= eval_data.r.const_trans;
                  break;
                case -1:
                  eval_data.r.const_trans = __data.coeffs[0]*__rval.const_trans;
                  eval_data.r.const_trans = 1./eval_data.r.const_trans;
                  break;
                case -2:
                  eval_data.r.const_trans = __data.coeffs[0]*__rval.const_trans;
                  eval_data.r.const_trans =
                            1./(eval_data.r.const_trans*eval_data.r.const_trans);
                  break;
                default:
                  eval_data.r.const_trans =
                                  pow(__data.coeffs[0]*__rval.const_trans,
                                      __data.params.nn());
                  break;
              }
            }
            else if(__data.params.nn() == 0)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = 1;
            }
            else if(__data.params.nn() != 1)
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_POW:
            if(eval_data.n == 0)
            {
              eval_data.r.was_lin = __rval.was_lin;
              eval_data.r.node_num = __rval.node_num;
            }
            if(__rval.was_const)
            {
              if(eval_data.n == 0)
              {
                eval_data.r.const_trans = __rval.const_trans;
                eval_data.r.was_const = true;
              }
              else if(eval_data.n == 1)
              {
                if(eval_data.r.was_const)
                {
                  eval_data.r.const_trans = pow(__data.coeffs[0]*
                      eval_data.r.const_trans+__data.params.nd(),
                      __data.coeffs[1]*__rval.const_trans);
                  eval_data.r.was_lin = true;
                }
                else if(__rval.const_trans == 0)
                {
                  eval_data.r.const_trans = 1;
                  eval_data.r.was_lin = eval_data.r.was_const = true;
                }
                else if(__rval.const_trans == 1 && eval_data.r.was_lin
                        && eval_data.old_nlin == -1)
                { // this is bad, since that means that we have to
                  // rewalk the left part of the graph
                  lincoeff_visitor lv(eval_data.ranges, eval_data.coeffs,
                                      eval_data.constant, eval_data.mod);
                  lv.eval_data.trans = eval_data.old_trans;
                  eval_data.r = evaluate(lv,
                                  eval_data.mod->node(eval_data.r.node_num));
                }
                else
                  eval_data.r.was_lin = eval_data.r.was_const = false;
              }
            }
            else
              eval_data.r.was_const = eval_data.r.was_lin = false;
            if(eval_data.n == 1)
              eval_data.r.node_num = __data.node_num;
            break;
          case EXPRINFO_EXP:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = exp(__data.coeffs[0]*__rval.const_trans+
                  __data.params.nd());
            }
            else
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_LOG:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = log(__data.coeffs[0]*__rval.const_trans+
                  __data.params.nd());
            }
            else
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_SIN:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = sin(__data.coeffs[0]*__rval.const_trans+
                  __data.params.nd());
            }
            else
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_COS:
            if(__rval.was_const)
            {
              eval_data.r.was_const = eval_data.r.was_lin = true;
              eval_data.r.const_trans = cos(__data.coeffs[0]*__rval.const_trans+
                  __data.params.nd());
            }
            else
              eval_data.r.was_lin = false;
            break;
          case EXPRINFO_ATAN2:
            if(__rval.was_const)
            {
              if(eval_data.n == 0)
              {
                eval_data.r.const_trans = __rval.const_trans;
                eval_data.r.was_const = true;
              }
              else if(eval_data.n == 1)
              {
                if(eval_data.r.was_const)
                {
                  eval_data.r.const_trans =
                          atan2(eval_data.r.const_trans, __rval.const_trans);
                  eval_data.r.was_lin = true;
                }
                else
                  eval_data.r.was_lin = eval_data.r.was_const = false;
              }
            }
            else
              eval_data.r.was_const = eval_data.r.was_lin = false;
            break;
        case EXPRINFO_GAUSS:
          if(__rval.was_const)
          {
            eval_data.r.was_const = eval_data.r.was_lin = true;
            eval_data.r.const_trans =
              (__data.coeffs[0]*__rval.const_trans-__data.params.d()[0])/
                                              __data.params.d()[1];
            eval_data.r.const_trans = exp(eval_data.r.const_trans*
                eval_data.r.const_trans);
          }
          else
            eval_data.r.was_lin = false;
          break;
        case EXPRINFO_LIN:
        case EXPRINFO_QUAD:
          std::cerr << "Must not be here in lincoeff_visitor" << std::endl;
          throw "Strange error";
          break;
        default:
          std::cerr << "Unknown expression type in get_linear_coeffs" <<
            std::endl;
          throw "Programming error";
          break;
      }
    }
    else if(__data.operator_type > 0)
    {
      eval_data.r.was_lin = false;
    }
    eval_data.r.was_lin = eval_data.r.was_lin || __rval.was_lin;
    if(++eval_data.n < __data.n_children)
      eval_data.trans = eval_data.old_trans * __data.coeffs[eval_data.n];
    return 0;
  }
  
  int update(const lincoeff_visitor_ret& __rval)
        {
          if(__rval.was_const)
            *eval_data.constant += __rval.const_trans;
          eval_data.r.was_lin = eval_data.r.was_lin || __rval.was_lin;
          return 0;
        }

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

inline bool model::get_linear_coeffs(const walker& expr, sparse_vector<double>& coeffs,
                        double& constant, const std::vector<interval>& _ranges)
{
  lincoeff_visitor lv(&_ranges, &coeffs, &constant, this);
  lincoeff_visitor_ret r = evaluate(lv, expr);
  if(r.was_const)
    constant += r.const_trans;
  return r.was_lin;
}

inline bool model::get_linear_coeffs(const walker& expr, sparse_vector<double>& coeffs,
                              double& constant)
{
  lincoeff_visitor lv(NULL, &coeffs, &constant, this);
  lincoeff_visitor_ret r = evaluate(lv, expr);
  if(r.was_const)
    constant += r.const_trans;
  return r.was_lin;
}

inline void model::simplify_visitor_0::postorder_help(const expression_node &r,
    unsigned int n_chld)
{
  if(r.operator_type <= 0)
  {
    switch(r.operator_type)
    {
      case EXPRINFO_GHOST:
        { model::walker _w = __mod->node(r.node_num);
          postorder_help(*_w, _w.n_children());     // this sets flags, m,...
          v_i = _w->var_indicator();
          s = _w->sem;
        }
        flags |= SIMPLIFY_0_IS_GHOST;
        break;
      case EXPRINFO_CONSTANT:
        s.degree = 0;
        s.dim = 0;
        s.property_flags.c_info = c_linear;
        s.property_flags.separable = t_true;
        flags |= SIMPLIFY_0_IS_CONST;
        if(is_integer(r.params.nd()))
        {
          int i = get_integer(r.params.nd());
          flags |= SIMPLIFY_0_CONST_IS_INTEGER;
          m = i;
        }
        else
          m = r.params.nd();
        break;
      case EXPRINFO_VARIABLE:
        s.degree = 1;
        s.dim = 1;
        s.property_flags.c_info = c_linear;
        s.property_flags.separable = t_true;
        flags |= SIMPLIFY_0_IS_VAR;
        v_i.set(r.params.nn());
        break;
      case EXPRINFO_SUM:
        flags |= SIMPLIFY_0_IS_SUM;
        m = r.params.nd();
        transfer = &r.coeffs;
        if(n_chld == 1) // this is \la+<node>
          flags |= SIMPLIFY_0_SUM_IS_SIMPLE;
        break;
      case EXPRINFO_PROD:
        flags |= SIMPLIFY_0_IS_PROD;
        m = r.params.nd();
        transfer = &r.coeffs;
        if(n_chld == 1) // this is \la*<node>
          flags |= SIMPLIFY_0_PROD_IS_SIMPLE;
        break;
      case EXPRINFO_SQUARE:
      case EXPRINFO_INTPOWER:
        break;
      case EXPRINFO_INVERT:
      case EXPRINFO_DIV:
        flags |= SIMPLIFY_0_IS_MULTIPLICATIVE;
        m = r.params.nd();
        s.degree = -1;
        break;
      case EXPRINFO_IN:
      case EXPRINFO_IF:
      case EXPRINFO_AND:
      case EXPRINFO_OR:
      case EXPRINFO_NOT:
      case EXPRINFO_IMPLIES:
      case EXPRINFO_COUNT:
      case EXPRINFO_ALLDIFF:
      case EXPRINFO_LEVEL:
      case EXPRINFO_NEIGHBOR:
      case EXPRINFO_NOGOOD:
        s.annotation_flags.integer = true;
        // no break!
      case EXPRINFO_MAX:
      case EXPRINFO_MIN:
      case EXPRINFO_SQROOT:
      case EXPRINFO_ABS:
      case EXPRINFO_POW:
      case EXPRINFO_EXP:
      case EXPRINFO_LOG:
      case EXPRINFO_SIN:
      case EXPRINFO_COS:
      case EXPRINFO_ATAN2:
      case EXPRINFO_GAUSS:
        s.degree = -1;
        break;
      case EXPRINFO_LIN:
        break;
      case EXPRINFO_QUAD:
        break;
      default:
        // give bad messages
        break;
    }
  }
  else if(r.operator_type > 0)
  {
    s.degree = -1;
    s.property_flags.c_info = c_maybe;
    s.property_flags.separable = t_maybe;
  }
  s.dim = v_i.sum(0, num_of_vars);
  // now adapt semantics
  nnum = r.node_num;
  if(s.degree < 0) s.degree = -1;
}

//#include <model.cc>

#endif /* _MODEL_INLINE_H_ */

---