model.hpp

Go to the documentation of this file.

// Model implementation -*- C++ -*-

// $Id: model.hpp 875 2009-02-13 13:10:41Z dferi $
// 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_HPP_
#define _MODEL_HPP_

#include <api_exception.h>

namespace coco {

#define delnodes_insert(N) do {\
  std::vector<unsigned int>::iterator __x = \
    lower_bound((*delnodes).begin(), (*delnodes).end(), (unsigned int)(N)); \
    if(__x == (*delnodes).end() || *__x != (unsigned int)(N)) \
      (*delnodes).insert(__x, (unsigned int)(N)); \
  } while(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.empty())
  {
    std::pop_heap(node_free.begin(), node_free.end());
    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);
    std::push_heap(node_free.begin(), node_free.end());
  }
  while(!node_free.empty() && node_free[0] == node_num_max-1)
  {
    --node_num_max;
    std::pop_heap(node_free.begin(), node_free.end());
    node_free.pop_back();
  }
  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.empty())
  {
    std::pop_heap(var_free.begin(), var_free.end());
    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);
    std::push_heap(var_free.begin(), var_free.end());
  }
  while(!var_free.empty() && var_free[0] == var_num_max-1)
  {
    --var_num_max;
    std::pop_heap(var_free.begin(), var_free.end());
    var_free.pop_back();
  }
  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.empty())
  {
    std::pop_heap(const_free.begin(), const_free.end());
    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);
    std::push_heap(const_free.begin(), const_free.end());
  }
  while(!const_free.empty() &&
        const_free[0] == const_num_max-1)
  {
    --const_num_max;
    std::pop_heap(const_free.begin(), const_free.end());
    const_free.pop_back();
  }
  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::model_name() const
{
  if(mod_name.empty())
    return std::string("%%model");
  else
    return mod_name;
}

inline void model_iddata::model_name(const std::string& vn)
{
  mod_name = vn;
}

inline const std::string model_iddata::var_name(unsigned int n) const
{
  if(n >= var_names.size() || var_names[n].empty())
  {
    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].empty())
  {
    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.empty())
    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())
  {
    throw api_exception(apiee_internal, "index error in fixed_var!");
  }
  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())
  {
    throw api_exception(apiee_internal, "index error in unused_var!");
  }
  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())
  {
    throw api_exception(apiee_internal, "index error in unused_constr!");
  }
  return unused_constrs[n];
}

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

inline model_gid::model_gid(model& __m, model_iddata* __i)
        : model_ref(__m), glob_ref(), gvar_ref(), gconst_ref(), const_back_ref()
{
  if(__i == NULL)
    _id = new model_iddata;
  else
    _id = __i;
  _id->new_ref(*this);
}

inline model_gid::model_gid(model& __m, unsigned int n, model_iddata* __i)
         : model_ref(__m), glob_ref(), gvar_ref(), gconst_ref(), const_back_ref()
{
  if(__i == NULL)
    _id = new model_iddata(n);
  else
    _id = __i;
  _id->new_ref(*this);
}

inline model_gid::model_gid(model& __mr, const model_gid& __m)
                : _id(__m._id), model_ref(__mr), glob_ref(__m.glob_ref),
                  gvar_ref(__m.gvar_ref), gconst_ref(__m.gconst_ref),
                  const_back_ref(__m.const_back_ref)
{
  _id->new_ref(*this);
  std::vector<model::walker>::iterator __x, __e;
  model::walker mgr(__m.model_ref.ground()), gr(__mr.ground());
  __e = glob_ref.end();
  for(__x = glob_ref.begin(); __x != __e; ++__x)
    if(*__x == mgr)
      *__x = gr;
  __e = gvar_ref.end();
  for(__x = gvar_ref.begin(); __x != __e; ++__x)
    if(*__x == mgr)
      *__x = gr;
  __e = gconst_ref.end();
  for(__x = gconst_ref.begin(); __x != __e; ++__x)
    if(*__x == mgr)
      *__x = gr;
}

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 = (unsigned int)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)
{
  unsigned 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 api_exception(apiee_internal,
          "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 && gconst_ref[_n] != model_ref.ground())
  {
    __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 const std::map<int, std::vector<double > >& model_gid::trace_points()
        const { return _id->trace_points; }

inline const std::map<int, std::vector<double > >& model_gid::optima()
        const { return _id->optima; }

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),
                                destroy_node_ref(false), lin(), matd(),
                                mati(), ocoeff(0), objective(ground()),
                                constraints(), has_read_errors(false)
{ gid = new model_gid(*this, __num_of_vars); }

inline model::model(model_gid* __id, bool clone) :
                                       _Base(), node_ref(), var_ref(),
                                       ghost_ref(), destroy_node_ref(false),
                                       lin(), matd(), mati(),
                                       ocoeff(0), objective(ground()),
                                       constraints(), has_read_errors(false)
{
  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),
                                 destroy_node_ref(false),
                                 lin(__m.lin), matd(__m.matd),
                                 mati(__m.mati), ocoeff(__m.ocoeff),
                                 constraints(),
                                 has_read_errors(__m.has_read_errors)
{
#if MODEL_INLINE_DEBUG
  std::cerr << "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), destroy_node_ref(false),
                                 lin(__m.lin), matd(__m.matd),
                                 mati(__m.mati), ocoeff(__m.ocoeff),
                                 constraints(),
                                 has_read_errors(__m.has_read_errors)
{
  if(!__id)
  {
    throw api_exception(apiee_internal,
        "NULL gid is forbidden in model constructor!");
  } 
  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;
#endif
  if(destroy_node_ref)
  {
#if MODEL_INLINE_DEBUG
    std::cerr << "\tremoving node references" << std::endl;
#endif
    free_all();
  }
  if(keep_gid)
  {
#if MODEL_INLINE_DEBUG
    std::cerr << "\tremoving gid" << std::endl;
#endif
    delete gid;
  }
}

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 void model::free_nodes_on_destroy()
        {
#if MODEL_INLINE_DEBUG
          std::cerr << "\t(before) destroy_node_ref = " << destroy_node_ref << std::endl;
#endif
          destroy_node_ref = true;
#if MODEL_INLINE_DEBUG
          std::cerr << "\t(after) destroy_node_ref = " << destroy_node_ref << std::endl;
#endif
        }

inline void model::keep_nodes_on_destroy()
        {
#if MODEL_INLINE_DEBUG
          std::cerr << "\t(before) destroy_node_ref = " << destroy_node_ref << std::endl;
#endif
          destroy_node_ref = false;
#if MODEL_INLINE_DEBUG
          std::cerr << "\t(after) destroy_node_ref = " << destroy_node_ref << std::endl;
#endif
        }

inline const std::string model::model_name() const { return gid->model_name(); }
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 const std::map<int, std::vector<double > >& model::trace_points() const
        { return gid->trace_points(); }
inline const std::map<int, std::vector<double > >& model::optima() const
        { return gid->optima(); }
                                            
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_CONST_IS_EQUATION    (1<<25)

#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)

#define SIMPLIFY_0_SKIP_THIS_NODE       (1<<30)

class model::simplify_visitor_0 : public postorder_visitor<expression_node,
                                                simplify_visitor_0,
                                                const simplify_visitor_0&>
{
  int num_of_vars;
  semantics s;
  int n, ndel, ncorr;
  variable_indicator v_i;
  std::vector<unsigned int>* delnodes;
  std::vector<unsigned int>* delghosts;
  std::vector<std::pair<unsigned int,unsigned int> >* deledges;
  std::vector<std::pair<unsigned int,unsigned int> > trans_ghosts;
  std::vector<std::triple<unsigned int,unsigned int,unsigned int> >
                                                                  trans_edges;
  std::vector<bool>* visited;
  unsigned int flags;
  basic_alltype 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), ncorr(0), v_i(0),
                        delnodes(NULL), delghosts(NULL), deledges(NULL),
                        trans_ghosts(), trans_edges(), visited(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<unsigned int>* __dg,
                     std::vector<std::pair<unsigned int,unsigned int> >* __de,
                     std::vector<bool>* __vt, model* __m) :
                        num_of_vars(__nv), s(), n(0), ndel(0), ncorr(0),
                        v_i(__nv), delnodes(__dn), delghosts(__dg),
                        deledges(__de), trans_ghosts(), trans_edges(),
                        visited(__vt), 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), ncorr(__x.ncorr), v_i(__x.v_i),
                      delnodes(__x.delnodes), delghosts(__x.delghosts),
                      deledges(__x.deledges), trans_ghosts(),
                      trans_edges(),
                      visited(__x.visited), 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; ncorr = 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);

  bool preorder(expression_node &r)
  { 
    // check whether this node is already 'virtually' deleted
    std::vector<unsigned int>::iterator __x;
    if(r.operator_type == EXPRINFO_GHOST)
    {
      __x = lower_bound((*delghosts).begin(), (*delghosts).end(),
        r.node_num);
      if(__x != (*delghosts).end() && *__x == r.node_num)
      {
        flags |= SIMPLIFY_0_SKIP_THIS_NODE;
        return false;
      }
    }
    else
    {
      __x = lower_bound((*delnodes).begin(), (*delnodes).end(),
        r.node_num);
      if(__x != (*delnodes).end() && *__x == r.node_num)
      {
        flags |= SIMPLIFY_0_SKIP_THIS_NODE;
        return false;
      }
    }

    init(); 
    if((*visited)[r.node_num])
    {
      n = r.n_children;
      ndel = n-r.coeffs.size();
      v_i.set_all(r.v_i);
      s = r.sem;
      return false;
    }
    else
      return true;
  }
  
  simplify_visitor_0 value() { return *this; }
  simplify_visitor_0 vvalue() { s.degree = -2; return *this; }

  void simple_sum_prod_update(expression_node& r,
                              const simplify_visitor_0& __s)
  {
    std::vector<unsigned int>::iterator __x;
    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)
        trans_ghosts.push_back(std::make_pair(__s.nnum, r.node_num));
      else
      {
        model::walker _node = __mod->node(__s.nnum);
        if(_node->n_parents == 1)
          delnodes_insert(__s.nnum);
        else
        {
          trans_edges.push_back(std::make_triple(__s.nnum, r.node_num, n));
          _node->n_parents--;
        }
      }
    }
    else 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.;
    }
    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)
        trans_ghosts.push_back(std::make_pair(__s.nnum, r.node_num));
      else
      {
        model::walker _node = __mod->node(__s.nnum);
        if(_node->n_parents == 1)
          delnodes_insert(__s.nnum);
        else
        {
          trans_edges.push_back(std::make_triple(__s.nnum, r.node_num, n));
          _node->n_parents--;
        }
      }
    }
  }

  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;
    std::vector<unsigned int>::iterator __x;

    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);
    if(_go.n_parents() == 1)
    {
      // collect disconnected ghosts (will be disconnected below)
      __x = lower_bound((*delghosts).begin(), (*delghosts).end(), nnum);
      if(__x == (*delghosts).end() || *__x != nnum)
        (*delghosts).insert(__x, nnum);
    }
    for(unsigned int i=0; i<unum.size(); ++i)
    {
      if(i == 0)
      {
        if(__mod->exchange_ghost(_pr, _go, unum[i]))
          __mod->remove_edge_and_deattach(_pr, _go);
      }
      else
        __mod->safe_ghost(_pr, unum[i]);
      __x = lower_bound((*delghosts).begin(), (*delghosts).end(), unum[i]);

      // remove reconnected ghosts
      if(__x != (*delghosts).end() && *__x == nnum)
        (*delghosts).erase(__x);
    }
  }

  void vcollect(const simplify_visitor_0& __s);

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

class model::simplify_visitor_m : public postorder_visitor<expression_node,
                                                const simplify_visitor_m&,
                                                const simplify_visitor_m&>
{
  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;
  std::vector<bool>* visited;
  unsigned int flags;
  basic_alltype 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_m() : num_of_vars(0), s(), n(0), ndel(0), v_i(0),
                        delnodes(NULL), deledges(NULL), visited(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_m(int __nv, std::vector<unsigned int>* __dn,
                     std::vector<std::pair<unsigned int,unsigned int> >* __de,
                     std::vector<bool>* __vt, model* __m) :
                        num_of_vars(__nv), s(), n(0), ndel(0), v_i(__nv),
                        delnodes(__dn), deledges(__de), visited(__vt),
                        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_m(const simplify_visitor_m& __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),
                      visited(__x.visited), 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);

  bool preorder(expression_node &r)
  { init(); 
    if((*visited)[r.node_num])
    {
      n = r.n_children;
      v_i = r.v_i;
      s = r.sem;
      return false;
    }
    else
      return true;
  }

  const simplify_visitor_m& value() { return *this; }
  const simplify_visitor_m& vvalue() { s.degree = -2; return *this; }

  void vcollect(const simplify_visitor_m& __s) {}

  void collect(expression_node &r, const simplify_visitor_m& __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 = COCO_INF; } 

  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 = COCO_INF;
          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 = std::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] = std::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 = std::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 = std::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 = std::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 = std::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 = std::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 = std::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 = std::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 = std::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_all()
{
#if MODEL_INLINE_DEBUG
  std::cerr << "Removing all var refs" << std::endl;
#endif
  std::vector<int> _all_vars;
  _all_vars.reserve(var_ref.size());

  for(std::vector<walker>::iterator __i = var_ref.begin();
      __i != var_ref.end(); ++__i)
    _all_vars.push_back((*__i)->params.nn());
  for(std::vector<int>::iterator __i = _all_vars.begin();
      __i != _all_vars.end(); ++__i)
    free_var_num(*__i);

#if MODEL_INLINE_DEBUG
  std::cerr << "Removing all const refs" << std::endl;
#endif
  std::vector<std::pair<walker, unsigned int> > _all_consts;
  _all_consts.reserve(constraints.size());

  for(std::vector<walker>::iterator __j = constraints.begin();
      __j != constraints.end(); ++__j)
  {
    unsigned int cnum;
    if(gid->get_const_num((*__j)->node_num, cnum))
      _all_consts.push_back(std::make_pair(*__j, cnum));
  }
  for(std::vector<std::pair<walker, unsigned int> >::iterator
                __k = _all_consts.begin();
      __k != _all_consts.end(); ++__k)
    free_const_num((*__k).first, (*__k).second);

#if MODEL_INLINE_DEBUG
  std::cerr << "Removing all node refs" << std::endl;
#endif
  std::vector<unsigned int> _all_nodes;
  _all_nodes.reserve(node_ref.size());

  for(std::vector<walker>::iterator __i = node_ref.begin();
      __i != node_ref.end(); ++__i)
    _all_nodes.push_back((*__i)->node_num);
  for(std::vector<unsigned int>::iterator __i = _all_nodes.begin();
      __i != _all_nodes.end(); ++__i)
    free_node_num(*__i);

}

inline void model::free_node_num(unsigned int _nnum)
{
#if MODEL_INLINE_DEBUG
  std::cerr << "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::free_var_num(int _vnum)
{
#if MODEL_INLINE_DEBUG
  std::cerr << "Removing var " << _vnum << std::endl;
#endif
  ref_iterator __x(lower_bound(var_ref.begin(),
                var_ref.end(), _vnum, __docompare_variables()));
  if(__x != var_ref.end() && (*__x)->params.nn() == _vnum)
    var_ref.erase(__x);
  gid->remove_var_id(_vnum);
}

inline void model::free_const_num(const walker& _w, unsigned int _cnum)
{
#if MODEL_INLINE_DEBUG
  std::cerr << "Removing constraint " << _cnum << std::endl;
#endif
  ref_iterator __x;
  for(__x = constraints.begin(); __x != constraints.end(); ++__x)
    if((*__x) == _w) break;
  if(__x != constraints.end())
    constraints.erase(__x);
  gid->remove_const_id(_cnum);
}

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(-COCO_INF,COCO_INF);
    __r = between_back(ground(), sky(), _en);
    ghost_ref.insert(__x, __r);
    return __r;
  }
}

inline bool model::exchange_ghost(walker prnt, walker gh_old,
                                  unsigned int _nnum)
{
  ref_iterator __x;
  bool ret(false);

  __x = lower_bound(ghost_ref.begin(), ghost_ref.end(), _nnum,
                      __docompare_nodes());
  if(__x != ghost_ref.end() && (*__x)->node_num == _nnum)
    replace_edge_to_child(prnt, gh_old, *__x);
  else // make a new ghost
  {
    walker __r;
    expression_node _en(EXPRINFO_GHOST, _nnum);
    _en.params = false;
    _en.f_bounds = interval(-COCO_INF,COCO_INF);

    children_iterator _cit;
    _cit = prnt._C_w_cur->get_childentry_iterator((void*)gh_old.node());
    __r = between(prnt, _cit, sky(), sky().parent_end(), _en);
    if(gh_old.n_parents() > 1)
      remove_edge(prnt, gh_old);
    else
      ret = true;
    ghost_ref.insert(__x, __r);
  }
  return ret;
}

inline void model::safe_ghost(walker prnt, 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)
  {
    add_edge_back(prnt, *__x);
  }
  else // make a new one
  {
    walker __r;
    expression_node _en(EXPRINFO_GHOST, _nnum);
    _en.params = false;
    _en.f_bounds = interval(-COCO_INF,COCO_INF);
    __r = between_back(prnt, sky(), _en);
    ghost_ref.insert(__x, __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, const basic_alltype& _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,
                           const basic_alltype& _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,
                          const basic_alltype& _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));
}

template <typename _TT>
inline void model::read_vector(std::vector<_TT>& hd, char *& cq,
                               const char *format, const char *message,
                               int ln)
{
  _TT help;

  cq += strspn(cq+1, " \t")+1;
  for(; cq; )
  {
    char *ct;

    ct = strsep(&cq, ",");
    ct += strspn(ct, " \t");
    if(sscanf(ct, format, &help) != 1)
    {
      throw api_exception(apiee_io, std::string("Malformed Input: Line ")+
                          convert_to_str(ln)+": "+message);
    }
    hd.push_back(help);
  }
}

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 = (unsigned int)_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 = (unsigned int)_x.n_parents();
  }
}

struct model::detect_0chain_visitor_st
{
  std::pair<unsigned int,unsigned int> r; 
  unsigned int *f;                        
  std::vector<unsigned int>* t;           
  std::vector<unsigned int>* d;           
  std::vector<unsigned int>* nn;          
  bool new_number;                        
};

class model::detect_0chain_visitor : public
  cached_forward_evaluator_base<model::detect_0chain_visitor_st,
                        expression_node, std::pair<unsigned int,unsigned int>,
                        model::const_walker>
{
private:
  typedef cached_forward_evaluator_base<model::detect_0chain_visitor_st,
                expression_node, std::pair<unsigned int,unsigned int>,
                model::const_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;  
  linalg::sparse_vector<double>* coeffs;     
  //mtl::compressed1D<double>* coeffs;
  double* constant;                  
  double old_trans;                     
  double trans;                         
  double mm_help;                       
  lincoeff_visitor_ret r;               
  int old_nlin;                         
  int in_nlin;                          
  unsigned int n;                       
};

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

public:
  lincoeff_visitor(const std::vector<interval>* _rg,
                   linalg::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);
  }

  virtual const_walker short_cut_to(const expression_node& __data)
  {
    return eval_data.mod->sky();
  }

  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:
          throw api_exception(apiee_internal, "Must not be here in lincoeff_visitor::calculate");
          break;
        default:
          throw api_exception(apiee_internal,
              std::string("Unknown expression type in lincoeff_visitor::calculate: ")+
              convert_to_str(__data.operator_type));
          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,
                      (const_walker) 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,
                     (const_walker) 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,
                      (const_walker) 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 = std::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 = std::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 =
                                std::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 = std::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,
                      (const_walker) 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 = std::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 = std::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 = std::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 = std::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 =
                        std::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 = std::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:
        throw api_exception(apiee_internal, "Must not be here in lincoeff_visitor::update");
        break;
      default:
        throw api_exception(apiee_internal,
            std::string("Unknown expression type in lincoeff_visitor::update: ")+
            convert_to_str(__data.operator_type));
        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 const_walker& expr,
                linalg::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 const_walker& expr,
                linalg::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);
          if(_w->operator_type != EXPRINFO_GHOST)
            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:
      case EXPRINFO_MEAN:
        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;
}

inline void model::simplify_visitor_m::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);
          if(_w->operator_type != EXPRINFO_GHOST)
            postorder_help(*_w, (unsigned int)_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:
      case EXPRINFO_PROD:
      case EXPRINFO_SQUARE:
      case EXPRINFO_INTPOWER:
      case EXPRINFO_INVERT:
      case EXPRINFO_DIV:
        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;
}

#if DEBUG_EXPRESSIONS
} // namespace coco
#include <g_algotest.h>
namespace coco {
class __test_expr_help
{
public:
  unsigned int operator()(const expression_node& __n) const
  {
    return __n.test();
  }
};

inline void consistency_test(const model& __m)
{
  recursive_consistency_test(__m.ground(), __test_expr_help());
}
#endif

} // namespace coco

#endif /* _MODEL_HPP_ */

---