---

model.cc

Go to the documentation of this file.

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

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

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

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

#include <model.h>

extern const char *expr_names[];

void model_iddata::compress_numbers(bool renum_vars,
                                               bool renum_consts)
{
  unsigned int il, iu, jl, ju, idx = node_num_max-1;
  std::vector<std::pair<unsigned int,unsigned int> > __v;
  
  if(node_free.size() != 0)
  {
    sort(node_free.begin(), node_free.end());
    __v.reserve(node_free.size());
    for(il = 0, iu = node_free.size()-1, jl = node_free[il], ju = node_free[iu];
        il <= iu;)
    {
      if(idx <= ju)
      {
        if(idx == ju)
          --idx;
        else
          node_num_max++; /* this is a correction for the subtraction later */
        if(iu == 0)
          break;
        --iu;
        ju = node_free[iu];
      }
      else
      {
        __v.push_back(std::make_pair(idx,jl));
        --idx;
        ++il;
        jl = node_free[il];
      }
    }
    node_num_max -= node_free.size();
    node_free.erase(node_free.begin(), node_free.end());
  }
  if(__v.size() != 0)
  {
    for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
        ++__b)
      (*__b)->renumber_nodes(__v);
  }
  __v.erase(__v.begin(), __v.end());

  if(renum_vars && var_free.size() != 0)
  {
    sort(var_free.begin(), var_free.end());
    __v.reserve(var_free.size());
    for(il = 0, iu = var_free.size()-1, jl = var_free[il], ju = var_free[iu];
        il <= iu;)
    {
      if(idx <= ju)
      {
        if(idx == ju)
          --idx;
        else
          var_num_max++; /* this is a correction for the subtraction later */
        if(iu == 0)
          break;
        --iu;
        ju = var_free[iu];
      }
      else
      {
        __v.push_back(std::make_pair(idx,jl));
        --idx;
        ++il;
        jl = var_free[il];
      }
    }
    var_num_max -= var_free.size();
    var_free.erase(var_free.begin(), var_free.end());
  }
  if(__v.size() != 0)
  {
    for(std::vector<std::pair<unsigned int,unsigned int> >::iterator __i = __v.begin();
        __i != __v.end(); ++__i)
      var_names[(*__i).second] = var_names[(*__i).first];
    var_names.erase(var_names.begin(), var_names.begin()+number_of_variables());
    for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
        ++__b)
      (*__b)->renumber_variables(__v);
  }
  __v.erase(__v.begin(), __v.end());

  if(renum_consts && const_free.size() != 0)
  {
    sort(const_free.begin(), const_free.end());
    __v.reserve(const_free.size());
    for(il = 0, iu = const_free.size()-1, jl = const_free[il],
        ju = const_free[iu]; il <= iu;)
    {
      if(idx <= ju)
      {
        if(idx == ju)
          --idx;
        else
          const_num_max++; /* this is a correction for the subtraction later */
        if(iu == 0)
          break;
        --iu;
        ju = const_free[iu];
      }
      else
      {
        __v.push_back(std::make_pair(idx,jl));
        --idx;
        ++il;
        jl = const_free[il];
      }
    }
    const_num_max -= const_free.size();
    const_free.erase(const_free.begin(), const_free.end());
  }
  if(__v.size() != 0)
  {
    for(std::vector<std::pair<unsigned int,unsigned int> >::iterator __i = __v.begin();
        __i != __v.end(); ++__i)
      const_names[(*__i).second] = const_names[(*__i).first];
    const_names.erase(const_names.begin(), const_names.begin()+number_of_constraints());
    for(std::list<model_gid*>::iterator __b = backref.begin(); __b != backref.end();
        ++__b)
      (*__b)->renumber_constraints(__v);
  }
}

void model_gid::renumber_nodes(
        std::vector<std::pair<unsigned int,unsigned int> >& __v)
{
  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;
    glob_ref[i] = glob_ref[(*__i).first];
    glob_ref[i]->node_num = i;
    __mf = const_back_ref.find((*__i).first);
    if(__mf != const_back_ref.end())
    {
      const_back_ref.insert(std::make_pair(i,(*__mf).second));
      const_back_ref.erase(__mf);
    }
  }
}

model::model(model_gid* __id, const erased_part& _ep, bool clone) :
                                       _Base(_ep), node_ref(), var_ref(),
                                       ghost_ref(), lin(), matd(), mati(),
                                       ocoeff(0), objective(ground()),
                                       constraints()
{
  if(__id)
  {
    if(clone)
    {
      gid = new model_gid(*this, *__id);
      keep_gid = true;
    }
    else
    {
      gid = __id;
      keep_gid = false;
    }
  }
  else
  {
    std::cerr << "Erased parts with NULL gid are forbidden in model constructor!" <<  std::endl;
    throw "Programming error";
  } 

  std::vector<walker> __n(number_of_nodes(), ground());
  std::vector<walker> __v(number_of_variables(), ground());
  std::vector<unsigned int> __c(number_of_constraints(), 0);
  ref_iterator __x;

  fill_arrays(ground(), __n, __v, ghost_ref);
  for(__x = __n.begin(); __x != __n.end(); ++__x)
    if(*__x != ground())
    {
      node_ref.push_back(*__x);
      gid->mk_globref((*__x)->node_num, *__x);
    }
  for(__x = __v.begin(); __x != __v.end(); ++__x)
    if(*__x != ground())
    {
      var_ref.push_back(*__x);
      gid->mk_gvarref((*__x)->params.nn(), *__x);
    }

  objective = gid->empty_reference();
  ocoeff = 0;
}

void model::fill_arrays(walker __w, std::vector<walker>& __n,
                        std::vector<walker>& __v, std::vector<walker>& __g)
{
  children_iterator __it, __e;
 
  if(__w.is_ground())
  {
    __e = __w.child_end();
    for(__it = __w.child_begin(); __it != __e; ++__it)
      fill_arrays(__w>>__it, __n, __v, __g);
  }
  else if(__w.is_sky())
    return;
  else
  {
    if(__n[__w->node_num] == ground())
    {
      if(__w->operator_type != EXPRINFO_GHOST)
        __n[__w->node_num] = __w;
      if(__w->operator_type == EXPRINFO_VARIABLE)
        __v[__w->params.nn()] = __w;
      else if(__w->operator_type == EXPRINFO_GHOST)
      {
        ref_iterator __x(lower_bound(__g.begin(), __g.end(), __w->node_num,
                         __docompare_nodes()));
        if(__x == __g.end() || *__x != __w)
          __g.insert(__x, __w);
      }
      __e = __w.child_end();
      for(__it = __w.child_begin(); __it != __e; ++__it)
        fill_arrays(__w>>__it, __n, __v, __g);
    }
  }
}

void model::copy_contents(model_gid* gid, const model& __m)
{
  std::vector<walker> __n(__m.number_of_nodes(), ground());
  std::vector<walker> __v(__m.number_of_variables(), ground());
  std::vector<int> __c(__m.number_of_constraints(), -1);
  ref_iterator __x;
  unsigned int i;

  for(unsigned int i = 0; i < __m.number_of_constraints(); ++i)
  {
    walker gc(__m.constraint(i));
    if(!__m.gid->empty(gc))
      __c[i] = (int) gc->node_num;
  }
  ghost_ref.reserve(__m.ghost_ref.size());
  fill_arrays(ground(), __n, __v, ghost_ref);
  node_ref.reserve(__m.number_of_managed_nodes());
  var_ref.reserve(__m.number_of_managed_variables());
  for(__x = __n.begin(); __x != __n.end(); ++__x)
    if(*__x != ground())
    {
      node_ref.push_back(*__x);
      gid->mk_globref((*__x)->node_num, *__x);
    }
  for(__x = __v.begin(); __x != __v.end(); ++__x)
    if(*__x != ground())
    {
      var_ref.push_back(*__x);
      gid->mk_gvarref((*__x)->params.nn(), *__x);
    }

  constraints.reserve(__m.constraints.size());
  for(const_ref_iterator __b = __m.constraints.begin();
      __b != __m.constraints.end(); ++__b)
    constraints.push_back(gid->node((*__b)->node_num));
  for(i = 0; i < __c.size(); ++i)
    if(__c[i] != -1)
      gid->mk_gconstref(i, gid->node((unsigned int)__c[i]));

  if(__m.is_empty(__m.objective))
    objective = gid->empty_reference();
  else
    objective = gid->node(__m.objective->node_num);
}

model::model(const char* name, bool do_simplify) : _Base(), node_ref(),
                                var_ref(), ghost_ref(), lin(), matd(), mati(),
                                ocoeff(0), objective(ground()), constraints()
{
#if MODEL_INLINE_DEBUG
  std::cout << "model_constructor: 1" << std::endl;
#endif
  std::ifstream inp(name);
#if MODEL_INLINE_DEBUG
  std::cout << "model_constructor: 2" << std::endl;
#endif
           
  if(!inp.is_open())
  {
    std::string str("Could not open file `");
    str += name;
    str += "': ";
    str += strerror(errno);
    throw (const char *)str.c_str();
  }
#if MODEL_INLINE_DEBUG
  std::cout << "model_constructor: 6" << std::endl;
#endif
  gid = new model_gid(*this);
  keep_gid = true;
#if MODEL_INLINE_DEBUG
  std::cout << "model_constructor: 7" << std::endl;
#endif
  read(inp);
#if MODEL_INLINE_DEBUG
  std::cout << "model_constructor: 8" << std::endl;
#endif
  if(do_simplify)
    basic_simplify();
  else
    set_counters();
#if MODEL_INLINE_DEBUG
  std::cout << "model_constructor: 9" << std::endl;
#endif
  prepare_after_read();
#if MODEL_INLINE_DEBUG
  std::cout << "model_constructor: 10" << std::endl;
#endif
}

#define throw_read(MSG) do { std::cerr << "Line " << ln << ": " << MSG;\
                             throw "Malformed Input!"; } while(0)
#define MT_DBL    1
#define MT_INTV   2
#define MT_INT    3
#define MT_STORED 8

interval model::read_interval(char *& cq, int ln)
{
  double lo, up;

  if(cq[0] != '[')
    throw_read("Interval specifier: no `['");
  else
  {
    char *cr;
    
    cq++;
    cr = strsep(&cq, ",");
    if(cq == NULL)
      throw_read("Interval specifier: no `,' inside interval");
    if(!strcmp(cr, "-I"))
      lo = -INFINITY;
    else if(!strcmp(cr, "I"))
      lo = INFINITY;
    else if(sscanf(cr, "%lf", &lo) != 1)
      throw_read("Interval specifier: not a wellformed interval lower bound");
    cr = strsep(&cq, "]");
    if(cq == NULL)
      throw_read("Interval specifier: missing closing `]'");
    if(!strcmp(cr, "-I"))
      up = -INFINITY;
    else if(!strcmp(cr, "I"))
      up = INFINITY;
    else if(sscanf(cr, "%lf", &up) != 1)
      throw_read("Interval specifier: not a wellformed interval upper bound");
    cr = strsep(&cq, ",");
  }
  return interval(lo,up);
}

void model::read(std::istream &i)
{
  char buf[1000], *cp, *cq;
  double coef;
  int ln = 0, do_mng, found_mng, found_varidx, found_sem, _hlp;
  unsigned int nnum1, nnum2;
  std::vector<walker> nodes;
  std::vector<unsigned int> v_var_idx;
  std::map<unsigned int, std::pair<int,void *> > matrix_def;
  interval hi;

  objective = ground();
  ocoeff = 0;
  try
  {
    do
    {
      do_mng = 1;
      found_mng = 0;
      found_sem = 0;
      found_varidx = 0;
      hi.set(-INFINITY, INFINITY);

      i.getline(buf, 999);
      ln++;
      cp = buf;

      if(!strncmp(cp, "<N>", 3)) // read the global information
      { 
        char c[3] = "xx";

        cp += 3;
        cp += strspn(cp, " \t");
        switch(cp[0])
        {
          case 'c':
            {
              cp += strspn(cp+1, " \t")+1;

              if(sscanf(cp, "%u %u", &nnum1, &nnum2) != 2)
                throw_read("Malformed Global Info: Constraint Number");
              if(nnum1 >= number_of_constraints())
                gid->number_of_constraints(nnum1+1);
              gid->make_const_back_ref(nnum2, nnum1);
            }
            break;
          case 'C':
            {
              int stl;
              cp += strspn(cp+1, " \t")+1;

              char *__hlp = new char[strlen(cp)+1];
              if(sscanf(cp, "%u '%s", &nnum1, __hlp) != 2 ||
                 __hlp[stl = strlen(__hlp)-1] != '\'')
                throw_read("Malformed Global Info: Constraint Name");
              __hlp[stl] = '\0';
              if(nnum1 >= number_of_constraints())
                gid->number_of_constraints(nnum1+1);
              gid->const_name(nnum1, __hlp);
              delete [] __hlp;
            }
            break;
          case 'F':
            {
              int stl;
              cp += strspn(cp+1, " \t")+1;

              char *__hlp = new char[strlen(cp)+1];
              if(sscanf(cp, "'%s %lf", __hlp, &coef) != 2 ||
                 __hlp[stl = strlen(__hlp)-1] != '\'')
                throw_read("Malformed Global Info: Variable Name");
              __hlp[stl] = '\0';
              gid->fixed_var(__hlp, coef);
              delete [] __hlp;
            }
            break;
          case 'M':
            cp += strspn(cp+1, " \t")+1;
            if(sscanf(cp, "%u m%c%c", &nnum1, c, c+1) != 3)
              throw_read("Malformed Global Info: Objective function");
            if(!strcmp("in", c))
              ocoeff = 1;
            else if(!strcmp("ax", c))
              ocoeff = -1;
            else if(!strcmp("00", c))
              ocoeff = 0;
            else
              throw_read("Malformed Global Info: Objective function is neither max nor min");
            if(ocoeff != 0)
              objective = nodes[nnum1];
            else
              objective = ground();
            break;
          case 'N':
            {
              int stl;
              cp += strspn(cp+1, " \t")+1;

              char *__hlp = new char[strlen(cp)+1];
              if(sscanf(cp, "%u '%s", &nnum1, __hlp) != 2 ||
                 __hlp[stl = strlen(__hlp)-1] != '\'')
                throw_read("Malformed Global Info: Variable Name");
              __hlp[stl] = '\0';
              gid->var_name(nnum1, __hlp);
              delete [] __hlp;
            }
            break;
          case 'O':
            {
              int stl;

              cp += strspn(cp+1, " \t")+1;
              cq = strsep(&cp, ":");
              char *__hlp = new char[strlen(cq)+1];

              if(sscanf(cq, "'%s", __hlp) != 1 ||
                 __hlp[stl = strlen(__hlp)-1] != '\'')
                throw_read("Malformed Global Info: Objective Info");
              else
              {
                __hlp[stl] = '\0';
                gid->obj_name(__hlp);
              }
              if(cp && *cp != '\0')
              {
                cp += strspn(cp, " \t");
                cq = strsep(&cp, ",");

                if(sscanf(cq, "d %lf", &coef) != 1)
                  throw_read("Malformed Global Info: Objective Info");
                else
                  gid->obj_adj(coef);
                if(cp && *cp != '\0')
                {
                  cp += strspn(cp, " \t");
                  if(sscanf(cp, "%lf", &coef) != 1)
                    throw_read("Malformed Global Info: Objective Info");
                  else
                    gid->obj_mult(coef);
                }
              }
              delete [] __hlp;
            }
            break;
          case 'V':
            cp += strspn(cp+1, " \t")+1;
            if(sscanf(cp, "%u", &nnum1) != 1)
              throw_read("Malformed Global Info: Num of Vars");
            gid->number_of_variables(nnum1);
#if MODEL_INLINE_DEBUG
            std::cout << "gid : number of variables = " <<
                         gid->number_of_variables() << std::endl;
#endif
            break;
          case 'S':
            {
              int stl;
              cp += strspn(cp+1, " \t")+1;
              char *__hlp = new char[strlen(cp)+1];
              if(sscanf(cp, "'%s", __hlp) != 1 ||
                 __hlp[stl = strlen(__hlp)-1] != '\'')
                throw_read("Malformed Global Info: Automatic Constraint");
              else
              {
                __hlp[stl] = '\0';
                gid->unused_constr(__hlp);
              }
              delete [] __hlp;
            }
            break;
          case 'U':
            {
              int stl;
              cp += strspn(cp+1, " \t")+1;
              char *__hlp = new char[strlen(cp)+1];
              if(sscanf(cp, "'%s", __hlp) != 1 ||
                 __hlp[stl = strlen(__hlp)-1] != '\'')
                throw_read("Malformed Global Info: Unused Variable");
              else
              {
                __hlp[stl] = '\0';
                gid->unused_var(__hlp);
              }
              delete [] __hlp;
            }
            break;
          default: 
            throw_read("Malformed Global Info: Wrong type");
            break;
        }
      }
      else if(!strncmp(cp, "<M>", 3)) // read a matrix definition
      {
        bool is_for_lin = false;
        unsigned int mat_num = 0;
        char dens, shape;
        
        cp+=3;
        cp+=strspn(cp, " \t");
        if(sscanf(cp, "L %u %u", &nnum1, &nnum2) == 2)
        {
          is_for_lin = true;
          dens = 'D';
          shape = 'R';
        }
        else if(sscanf(cp, "%u %c %c %u %u", &mat_num, &dens, &shape, &nnum1,
                                        &nnum2) != 5)
          throw_read("Malformed Matrix Definition");
        if(dens != 'D' && dens != 'I' && dens != 'N')
          throw_read("Malformed Matrix Definition: must be double, interval, or integer");
        if(shape != 'R' && shape != 'S' && shape != 'T')
          throw_read("Malformed Matrix Definition: must be rectangular, symmetric, or triangular");
        if(shape == 'S' || shape == 'T')
          throw_read("Matrix shapes not yet implemented!");
        if(is_for_lin)
          lin = matrix<double>(nnum1,nnum2);
        else
        {
          void *matrix_storage;
          int matrix_type;

          matrix_type = (dens == 'D' ? MT_DBL :
                                        (dens == 'I' ? MT_INTV : MT_INT));
          switch(dens)
          {
            case 'D':
              matrix_storage = (void *) new matrix<double>(nnum1,nnum2);
              break;
            case 'I':
              matrix_storage = (void *) new matrix<interval>(nnum1,nnum2);
              break;
            case 'N':
              matrix_storage = (void *) new matrix<int>(nnum1,nnum2);
              break;
          }
          matrix_def.insert(std::make_pair(mat_num, 
                                std::make_pair(matrix_type, matrix_storage)));
        }
      }
      else if(!strncasecmp(cp, "<R>", 3)) // read a matrix row
      {
        unsigned int mat_num = 0;
        bool dense = (cp[1] == 'R');
        std::vector<unsigned int> s_idx;
        unsigned int dim_r, dim_c;
        unsigned int idx, ridx;
        int matrix_type;
        void *matrix_storage;

        cp+=3+strspn(cp+3, " \t");
        cq = strchr(cp, ':');
        if(!cq)
          throw_read("Malformed Matrix Row Definition");
        *cq++ = '\0';
        if(cp[0] == 'L' && isspace(cp[1]))
        {
          cp += 1+strspn(cp+1, " \t");
          if(sscanf(cp, "%u", &nnum1) != 1)
            throw_read("Malformed Matrix Row Definition: row number?");
          dim_r = lin.nrows();
          dim_c = lin.ncols();
          matrix_storage = (void *)(&lin);
          matrix_type = MT_DBL;
        }
        else if(sscanf(cp, "%u %u", &mat_num, &nnum1) == 2)
        {
          std::pair<int,void*> x;
          std::map<unsigned int,std::pair<int,void*> >::iterator _x;

          _x = matrix_def.find(mat_num);
          if(_x == matrix_def.end())
            throw_read("Malformed Matrix Row Definition: matrix not yet defined!");
          x = (*_x).second;
          matrix_type = x.first;
          matrix_storage = x.second;
          if(matrix_type & MT_STORED)
            throw_read("Malplaced Matrix Row Definition: matrix has been used!");
          switch(matrix_type)
          {
            case MT_DBL:
              dim_r = ((matrix<double>*)matrix_storage)->nrows();
              dim_c = ((matrix<double>*)matrix_storage)->ncols();
              break;
            case MT_INTV:
              dim_r = ((matrix<interval>*)matrix_storage)->nrows();
              dim_c = ((matrix<interval>*)matrix_storage)->ncols();
              break;
            case MT_INT:
              dim_r = ((matrix<int>*)matrix_storage)->nrows();
              dim_c = ((matrix<int>*)matrix_storage)->ncols();
              break;
          }
        }
        else
          throw_read("Malformed Matrix Row Definition: matrix & row number?");
        if(nnum1 >= dim_r)
          throw_read("Malformed Matrix Row Definition: row number out of range!");
        cp = cq;
        if(!dense)
        {
          cp = strchr(cp, ':');
          if(!cp)
            throw_read("Malformed Matrix Row Definition: sparse needs two vectors");
          *cp++ = '\0';
          cq += strspn(cq, " \t");
          for(; cq; )
          {
            char *ct;

            ct = strsep(&cq, ",");
            ct += strspn(ct, " \t");
            if(sscanf(ct, "%u", &nnum2) != 1)
              throw_read("Malformed Matrix Row Definition: index vector definition");
            if(nnum2 >= dim_c)
              throw_read("Malformed Matrix Row Definition: column number out of range!");
            s_idx.push_back(nnum2);
          }
        }
        cp += strspn(cp, " \t");
        idx = 0;
        for(; cp; )
        {
          char *ct;

          if(dense)
          {
            if(idx >= dim_c)
              throw_read("Malformed Matrix Row Definition: value vector too long!");
            ridx = idx;
          }
          else
          {
            if(idx >= s_idx.size())
              throw_read("Malformed Matrix Row Definition: value vector longer than index vector!");
            ridx = s_idx[idx];
          }
          ct = strsep(&cp, ",");
          ct += strspn(ct, " \t");
          switch(matrix_type)
          {
            case MT_DBL:
              if(sscanf(ct, "%lf", &coef) != 1)
                throw_read("Malformed Matrix Row Definition: value vector definition - no integer!");
              (*(matrix<double>*)matrix_storage)(nnum1,ridx) = coef;
              break;
            case MT_INTV:
              hi = read_interval(ct, ln);
              (*(matrix<interval>*)matrix_storage)(nnum1,ridx) = hi;
              break;
            case MT_INT:
              if(sscanf(ct, "%u", &nnum2) != 1)
                throw_read("Malformed Matrix Row Definition: value vector definition - no integer!");
              (*(matrix<int>*)matrix_storage)(nnum1,ridx) = nnum2;
              break;
          }
          idx++;
        }
        if(dense && idx < dim_c)
          throw_read("Malformed Matrix Row Definition: value vector too short!");
        else if(!dense && idx < s_idx.size())
          throw_read("Malformed Matrix Row Definition: value vector shorter than index vector!");
      }
      else if(!strncmp(cp, "<H>", 3)) // read helper information
      {
        std::cerr << "Helper information is not yet implemented!" << std::endl;
        throw "NYI";
      }
      else if((_hlp = sscanf(cp, "<E> %u %u %lf", &nnum1, &nnum2, &coef)) == 3
              || _hlp == 2) // read an edge
      {
        add_edge_back(nodes[nnum1], nodes[nnum2]);
        remove_edge(nodes[nnum1], sky());
        remove_edge(ground(), nodes[nnum2]);
        nodes[nnum1]->coeffs.push_back(_hlp == 2 ? 1.0 : coef);
      }
      else if(sscanf(cp, "<%u>", &nnum1) == 1) // read a node definition
      {
#if MODEL_INLINE_DEBUG
        std::cout << "Read a node : nnum = " << nnum1 << std::endl;
#endif
        additional_info_u mng(0.0);
        semantics _sem;
        char _h;

        if(nnum1 >= nodes.size())
          nodes.insert(nodes.end(), nnum1-nodes.size()+1, ground());
        for(cp = strchr(cp, '>')+1; cp; )
        {
          cq = strsep(&cp, ":");
          cq += strspn(cq, " \t");
          if(cq[1] != ' ')
            _h = ' ';
          else
            _h = cq[0];
          switch(_h)
          {
            case 'V':
              if(sscanf(cq, "V %u", &nnum2) != 1)
                throw_read("Variable definition");
              else
              {
                if(nnum2 >= gid->number_of_variables())
                  new_variables(nnum2+1);
                
                expression_node _en(EXPRINFO_VARIABLE, nnum1);

                _en.v_i.reserve(gid->number_of_variables());
                _en.v_i.set(nnum2);
                _en.params = (int)nnum2;
                nodes[nnum1] = store_node(between_back(ground(), sky(), _en));

                gid->mk_gvarref(nnum2, nodes[nnum1]);
                do_mng = 0;
              }
              break;
            case 'C':
              if(sscanf(cq, "C %lf", &coef) != 1)
                throw_read("Constant definition");
              else
              {
                expression_node _en(EXPRINFO_CONSTANT, nnum1);
                _en.v_i.reserve(gid->number_of_variables());
                _en.v_i.clear();
                _en.params = coef;

                nodes[nnum1] = store_node(between_back(ground(), sky(), _en));
                do_mng = 0;
              }
              break;
            case 'G':
              if(sscanf(cq, "G %c", &_h) != 1 || (_h != 'T' && _h != 'F'))
                throw_read("Ghost definition");
              else
              {
                expression_node _en(EXPRINFO_GHOST, nnum1);
                _en.v_i.reserve(gid->number_of_variables());
                _en.v_i.clear();
                _en.params = (_h == 'T' ? true : false);
                _en.f_bounds = interval(-INFINITY,INFINITY);

                nodes[nnum1] = store_node(between_back(ground(), sky(), _en));
                do_mng = 0;
              }
              break;
            case 'l':
              {
                std::vector<bool> hd;
                char _l;

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

                  ct = strsep(&cq, ",");
                  ct += strspn(ct, " \t");
                  if(sscanf(ct, "%c", &_l) != 1 || (_l != 'T' && _l != 'F'))
                    throw_read("Boolean additional definition");
                  hd.push_back(_l == 'T' ? true : false);
                }
                if(hd.size() == 1)
                  mng = hd[0];
                else
                  mng = hd;
                found_mng = 1;
              }
              break;
            case 'S':
              {
                char *ct;

                cq += strspn(cq+1, " \t")+1;
                if(cq[0] != '\"')
                  throw_read("String additional definition");
                ct = cq+1;
                while(1)
                {
                  ct = strchr(ct, '\"');
                  if(ct == NULL)
                    throw_read("String additional definition");
                  if(*(ct-1) != '\\')
                    break;
                  ct += 1;
                }
                
                *ct = '\0';
                mng = cq;
                cq = ct+1;
                found_mng = 1;
              }
              break;
            case 'd':
              {
                std::vector<double> hd;

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

                  ct = strsep(&cq, ",");
                  ct += strspn(ct, " \t");
                  if(sscanf(ct, "%lf", &coef) != 1)
                    throw_read("Double additional definition");
                  hd.push_back(coef);
                }
                if(hd.size() == 1)
                  mng = hd[0];
                else
                  mng = hd;
                found_mng = 1;
              }
              break;
            case 'n':
              {
                std::vector<int> hd;
                int ihlp;

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

                  ct = strsep(&cq, ",");
                  ct += strspn(ct, " \t");
                  if(sscanf(ct, "%d", &ihlp) != 1)
                    throw_read("Integer additional definition");
                  hd.push_back(ihlp);
                }
                if(hd.size() == 1)
                  mng = hd[0];
                else
                  mng = hd;
                found_mng = 1;
              }
              break;
            case 'i':
              {
                std::vector<interval> hd;

                cq += strspn(cq+1, " \t")+1;
                for(; cq; )
                {
                  cq += strspn(cq, " \t");
                  hd.push_back(read_interval(cq, ln));
                }
                if(hd.size() == 1)
                  mng = hd[0];
                else
                  mng = hd;
                found_mng = 1;
              }
              break;
            case 'm':
              if(sscanf(cq, "m %u", &nnum2) != 1)
                throw_read("Matrix reference");
              else
              {
                do_mng = 2;
                found_mng = nnum2+2;
              }
              break;
            case 'v':
              {
                std::vector<unsigned int> hd;
                int ihlp;

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

                  ct = strsep(&cq, ",");
                  ct += strspn(ct, " \t");
                  if(sscanf(ct, "%d", &ihlp) != 1)
                    throw_read("Var index additional definition");
                  hd.push_back(ihlp);
                }
                found_varidx = 1;
                v_var_idx = hd;
              }
              break;
            case 'b':
              cq += strspn(cq+1, " \t")+1;
              hi = read_interval(cq, ln);
              break;
            case 's':
              cq += strspn(cq+1, " \t")+1;
              _sem.read(cq);
              found_sem = 1;
              break;
            default:
              for(nnum2 = 0; expr_names[nnum2]; nnum2++)
                if(!strcmp(expr_names[nnum2], cq))
                  break;
              if(expr_names[nnum2] == NULL)
                 throw_read("Unknown operator");
              else
              {
                expression_node _en(-nnum2, nnum1);
                
                _en.v_i.reserve(gid->number_of_variables());
                _en.v_i.clear();
#if MODEL_INLINE_DEBUG
                std::cout << "Expression node number = " << _en.node_num <<
                             std::endl;
#endif
                nodes[nnum1] = store_node(between_back(ground(), sky(), _en));
#if MODEL_INLINE_DEBUG
                std::cout << "Sky node number = " << sky()->node_num <<
                             std::endl;
#endif
              }
              break;
          }
        }
        if(!found_mng && (nodes[nnum1]->operator_type == EXPRINFO_PROD ||
           nodes[nnum1]->operator_type == EXPRINFO_INVERT))
        {
          mng = 1.;
          do_mng = 1;
        }
        if(do_mng == 2)
        {
          int matrix_type;
          void* matrix_storage;
          std::map<unsigned int,std::pair<int,void*> >::iterator __x;
          std::pair<int,void*> x;

          nnum2 = found_mng-2;
          __x = matrix_def.find(nnum2);
          if(__x == matrix_def.end())
            throw_read("Matrix reference is referencing undefined matrix");
          x = (*__x).second;
          matrix_type = x.first;
          matrix_storage = x.second;
          switch(matrix_type & ~MT_STORED)
          {
            case MT_DBL:
              if(matrix_type & MT_STORED)
                nodes[nnum1]->params = *(matrix<double>*)matrix_storage;
              else
                nodes[nnum1]->params.set_m((matrix<double>*)
                                                matrix_storage);
              break;
            case MT_INTV:
              if(matrix_type & MT_STORED)
                nodes[nnum1]->params = *(matrix<interval>*)matrix_storage;
              else
                nodes[nnum1]->params.set_im((matrix<interval>*)
                                                matrix_storage);
              break;
            case MT_INT:
              if(matrix_type & MT_STORED)
                nodes[nnum1]->params = *(matrix<int>*)matrix_storage;
              else
                nodes[nnum1]->params.set_nm((matrix<int>*)matrix_storage);
              break;
          }
          if(!(matrix_type & MT_STORED))
          {
            matrix_type |= MT_STORED;
            matrix_def.erase(__x);
            matrix_def.insert(std::make_pair(nnum2, std::make_pair(matrix_type,
                    matrix_storage)));
          }
        }
        else if(do_mng)
          nodes[nnum1]->params = mng;
        if(found_sem)
          nodes[nnum1]->sem = _sem;
        nodes[nnum1]->set_bounds(hi);
        nodes[nnum1]->n_parents = 1;
        if(found_varidx)
        {
          nodes[nnum1]->is_var = v_var_idx.size();
          nodes[nnum1]->var_idx = v_var_idx;
        }
      }
      else if(buf[0] != '\0')
        throw_read("Unknown input");
    } while(buf[0] != '\0');
  }
  catch(const char *str)
  {
    std::cerr << "Error: malformed input on line " << ln << ": " << str << "!" <<
         std::endl;
    std::terminate();
  }
  clr_sky_ground_link();
  for(std::map<unsigned int,std::pair<int,void*> >::iterator __x = matrix_def.begin();
      __x != matrix_def.end(); ++__x)
  {
    std::pair<int,void*> x;
    x = (*__x).second;
    if(!(x.first & MT_STORED))
    {
      switch(x.first)
      {
        case MT_DBL:
          delete (matrix<double>*) x.second;
          break;
        case MT_INTV:
          delete (matrix<interval>*) x.second;
          break;
        case MT_INT:
          delete (matrix<int>*) x.second;
          break;
      }
    }
  }
  for(int i = node_ref.size()-1; i >= 0; --i)
  {
    if(node_ref[i] == ground())
    {
      std::cerr << "Undefined node in node_ref!" << std::endl;
      throw "Programming error";
    }
    else
    {
      walker __w = node_ref[i] << node_ref[i].parent_begin();
      if(__w == ground())
        node_ref[i]->n_parents = 0;
    }
  }
#if MODEL_INLINE_DEBUG
  std::cout << "BrYcE : " << sky()->node_num << std::endl;
#endif
}

bool model::simplify_visitor_0::postorder(expression_node &r)
{
  postorder_help(r, n-ndel);
  // now adapt semantics
  if(r.sem.convexity().t() == 0)
    r.sem.convexity(s.convexity());
  r.sem.separability(s.separability());
  r.sem.degree = s.degree;
  r.sem.dim = s.dim;

  r.v_i = v_i;
  if(r_idx != -1)
  {
    int ne_del = n-ndel-r_idx-1;
    r.coeffs.erase(r.coeffs.begin()+r_idx+1,
                   r.coeffs.end());
    ndel += ne_del;
    r_idx = -1;
  }
  if(r.operator_type == EXPRINFO_CONSTANT ||
     r.operator_type == EXPRINFO_VARIABLE ||
     r.operator_type == EXPRINFO_GHOST)
    r.n_children = 0;
  else
    r.n_children = n - ndel;
  return false;
}
  
void model::simplify_visitor_0::vcollect(const simplify_visitor_0& __s)
{
  if((__s.flags & SIMPLIFY_0_IS_CONST)) // this is not a real constraint
  {
    double help = __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                           __s.m.nd();
    walker _n = __mod->node(__s.nnum);
    if(!_n->f_bounds.is_entire()) // this is not a pure objective
    {
      if(_n->f_bounds.contains(help)) // this constraint is always fulfilled
        (*delnodes).push_back(__s.nnum);
    }
    if(_n == __mod->objective) // the objective is constant
    {
      __mod->ocoeff = 0;
      __mod->gid->obj_adj(__mod->gid->obj_adj()+help);
      if(_n->operator_type == EXPRINFO_CONSTANT) // objective is really a const.
        _n->params = 0.;
    }
  }
  if((__s.flags & SIMPLIFY_0_IS_SUM &&
      __s.flags & SIMPLIFY_0_SUM_IS_SIMPLE) ||
     (__s.flags & SIMPLIFY_0_IS_PROD &&
      __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE))
  { // there is a l*<node>+m or l*<node> expression on top -> remove
    walker _n = __mod->node(__s.nnum);
    walker _c;
    std::vector<unsigned int>::const_iterator _b, _e((*delnodes).end());

    // find the only non-deleted child
    for(children_iterator _ci = _n.child_begin(); _ci != _n.child_end(); ++_ci)
    {
      _c = _n >> _ci;
      unsigned int _cnn = _c->node_num;
      for(_b = (*delnodes).begin(); _b != _e; ++_b)
        if((*_b) == _cnn)
          break;
      if(_b == _e) // this is it
        break;
    }
    
    if(!_n->f_bounds.is_entire())
    {
      if(__s.flags & SIMPLIFY_0_IS_SUM)
      {
        _n->f_bounds -= __s.m.nd();
        _n->f_bounds /= (*__s.transfer)[0];
      }
      else
        _n->f_bounds /= __s.m.nd();
      _c->f_bounds.intersect(_n->f_bounds);
    }
    if(_n == __mod->objective)
    {
      double a = __mod->gid->obj_adj(),
             m = __mod->gid->obj_mult();
      if(__s.flags & SIMPLIFY_0_IS_SUM)
      {
        a += m * __s.m.nd();
        m *= (*__s.transfer)[0];
      }
      else
        m *= __s.m.nd();
      __mod->gid->obj_adj(a);
      if(m < 0)
      {
        __mod->ocoeff = -__mod->ocoeff;
        __mod->gid->obj_adj(-__mod->obj_adj());
        m = -m;
      }
      __mod->gid->obj_mult(m);
    }
    (*delnodes).push_back(__s.nnum);
  }
  if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
  // we can transfer a factor
  {
    walker _n = __mod->node(__s.nnum);
    double help = __s.m.nd();
    _n->f_bounds /= help;
    if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
      _n->params = 1.;
    if(_n == __mod->objective)
    {
      if(help < 0)
      {
        __mod->ocoeff = -__mod->ocoeff;
        __mod->gid->obj_adj(-__mod->obj_adj());
        help = -help;
      }
      help *= __mod->gid->obj_mult();
      __mod->gid->obj_mult(help);
    }
  }
}

void model::simplify_visitor_0::collect(expression_node &r,
    const simplify_visitor_0& __s)
{
  nnum = r.node_num;
  if(r.operator_type < 0)
  {
    switch(r.operator_type)
    {
      case EXPRINFO_CONSTANT:
      case EXPRINFO_VARIABLE:
      case EXPRINFO_GHOST:
        // constants, variables, and ghosts have no children
        break;
      case EXPRINFO_SUM:
        if(n == 0)
        {
          s.degree = __s.s.degree;
          s.property_flags = __s.s.property_flags;
          if(r.coeffs[0] < 0)
            s.property_flags.c_info == -s.property_flags.c_info;
          if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]) &&
             is_integer(r.params.nd()))
            s.annotation_flags.integer = true;
          r_idx = -1;
        }
        else
        { 
          if(s.degree == -1 || __s.s.degree == -1)
            s.degree = -1;
          else if(__s.s.degree > s.degree)
            s.degree = __s.s.degree;
          s.property_flags.c_info = convex_e();
          if(s.property_flags.separable != t_true ||
             __s.s.property_flags.separable != t_true)
            s.property_flags.separable = t_maybe;
          if(s.annotation_flags.integer && 
             !(__s.s.annotation_flags.integer && is_integer(r.coeffs[n])))
            s.annotation_flags.integer = false;
          if(l_nnum == __s.nnum)
          {
            if(r_idx == -1)
              r_idx = n-ndel-1;
            r.coeffs[r_idx] += r.coeffs[n-ndel];
            (*deledges).push_back(std::make_pair(nnum, __s.nnum));
          }
          else if(r_idx != -1)
          {
            int ne_del = n-ndel-r_idx-1;
            r.coeffs.erase(r.coeffs.begin()+r_idx+1,
                           r.coeffs.begin()+n-ndel);
            ndel += ne_del;
            r_idx = -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() + 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();
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.coeffs[n-ndel] *= __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        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>::const_iterator b = __s.transfer->begin();
              b != __s.transfer->end(); ++b)
            r.coeffs.push_back(*b * help);
          n += __s.transfer->size();
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        break;
      case EXPRINFO_MEAN:
        if(n == 0)
        {
          s.degree = __s.s.degree;
          s.property_flags = __s.s.property_flags;
          if(r.coeffs[0] < 0)
          {
            r.operator_type = EXPRINFO_SUM;
            s.property_flags.c_info == -s.property_flags.c_info;
          }
          r_idx = -1;
        }
        else
        { 
          if(s.degree == -1 || __s.s.degree == -1)
            s.degree = -1;
          else if(__s.s.degree > s.degree)
            s.degree = __s.s.degree;
          s.property_flags.c_info = convex_e();
          if(s.property_flags.separable != t_true ||
             __s.s.property_flags.separable != t_true)
            s.property_flags.separable = t_maybe;
          if(l_nnum == __s.nnum)
          {
            if(r_idx == -1)
              r_idx = n-ndel-1;
            r.coeffs[r_idx] += r.coeffs[n-ndel];
            (*deledges).push_back(std::make_pair(nnum, __s.nnum));
          }
          else if(r_idx != -1)
          {
            int ne_del = n-ndel-r_idx-1;
            r.coeffs.erase(r.coeffs.begin()+r_idx+1,
                           r.coeffs.begin()+n-ndel);
            ndel += ne_del;
            r_idx = -1;
          }
        }
        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_SUM;
          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.operator_type = EXPRINFO_SUM;
          r.coeffs[n-ndel] *= __s.m.nd();
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.operator_type = EXPRINFO_SUM;
          r.coeffs[n-ndel] *= __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        else if(__s.flags & SIMPLIFY_0_IS_SUM) 
        {
          double help = r.coeffs[n-ndel];

          r.operator_type = EXPRINFO_SUM;
          r.params = r.params.nd() + __s.m.nd();
          r.coeffs.erase(r.coeffs.begin()+(n-ndel));
          ndel++;
          for(std::vector<double>::const_iterator b = __s.transfer->begin();
              b != __s.transfer->end(); ++b)
            r.coeffs.push_back(*b * help);
          n += __s.transfer->size();
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        break;
      case EXPRINFO_PROD:
        if(n == 0)
        {
          s.degree = __s.s.degree;
          s.property_flags = __s.s.property_flags;
          m = 1.;
          if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]) &&
             is_integer(r.params.nd()))
            s.annotation_flags.integer = true;
        }
        else
        { 
          if(__s.s.degree == -1 || s.degree == -1)
            s.degree = -1;
          else
            s.degree += __s.s.degree;
          s.property_flags.c_info = convex_e();
          s.property_flags.separable = t_maybe;
          if(s.annotation_flags.integer && 
             !(__s.s.annotation_flags.integer && is_integer(r.coeffs[n])))
            s.annotation_flags.integer = false;
        }
        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;
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.params = r.params.nd() * __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        if(r.params.nd() != 1)
        {
          flags |= SIMPLIFY_0_IS_MULTIPLICATIVE;
        }
        break;
      case EXPRINFO_MONOME:
        if(n == 0)
        {
          s.degree = __s.s.degree == -1 ? -1 : __s.s.degree*r.params.n()[0];
          s.property_flags = __s.s.property_flags;
          s.property_flags.c_info = convex_e();
          m = 1.;
          if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]))
            s.annotation_flags.integer = true;
        }
        else
        { 
          if(__s.s.degree == -1 || s.degree == -1)
            s.degree = -1;
          else
            s.degree += __s.s.degree*r.params.n()[n];
          s.property_flags.c_info = convex_e();
          s.property_flags.separable = t_maybe;
          if(s.annotation_flags.integer && 
             !(__s.s.annotation_flags.integer && is_integer(r.coeffs[n])))
            s.annotation_flags.integer = false;
        }
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { 
          double help =
              __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                           __s.m.nd();
          m = m.nd() * pow(help, r.params.n()[n]);
          r.coeffs.pop_back();
          (*delnodes).push_back(__s.nnum);
          ndel++;
        }
        else if(__s.flags & SIMPLIFY_0_IS_PROD)
        {
          m = m.nd()*pow(__s.m.nd(),r.params.n()[n]);
          r.coeffs.insert(r.coeffs.end(), __s.transfer->size()-1, 1.0);
          r.params.n().insert(r.params.n().begin()+n, __s.transfer->size()-1,
              r.params.n()[n]);
          n += __s.transfer->size()-1;
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          m = m.nd() * pow(__s.m.nd(),r.params.n()[n]);
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        if(m.nd() != 1)
        {
          flags |= (SIMPLIFY_0_IS_MULTIPLICATIVE | SIMPLIFY_0_IS_CORRECTEDMULT);
        }
        break;
      case EXPRINFO_MAX:
      case EXPRINFO_MIN:
        if(n == 0)
        {
          if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]) &&
             is_integer(r.params.nd()))
            s.annotation_flags.integer = true;
        }
        else if(s.annotation_flags.integer && 
             !(__s.s.annotation_flags.integer && is_integer(r.coeffs[n])))
          s.annotation_flags.integer = false;
        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.;
        }
        break;
      case EXPRINFO_SCPROD:
        if(n == 0)
        {
          s.property_flags.c_info = convex_e();
          if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]))
            s.annotation_flags.integer = true;
        }
        else if(n & 1) // n odd
        { 
          r.coeffs[n-1] *= r.coeffs[n];
          r.coeffs[n] = 1.;
          if(s.degree == -1 || __s.s.degree == -1 || l_s.degree == -1)
            s.degree = -1;
          else if(__s.s.degree+l_s.degree > s.degree)
            s.degree = __s.s.degree+l_s.degree;
          s.property_flags.separable = t_maybe;
          if(s.annotation_flags.integer &&
             !(__s.s.annotation_flags.integer && is_integer(r.coeffs[n-1])))
            s.annotation_flags.integer = false;
        }
        else
        {
          if(s.annotation_flags.integer &&
             !(__s.s.annotation_flags.integer && is_integer(r.coeffs[n])))
            s.annotation_flags.integer = false;
        }
        if(__s.flags & SIMPLIFY_0_IS_PROD &&
                __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE)
        {
          r.coeffs[n-ndel] *= __s.m.nd();
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.coeffs[n-ndel] *= __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        break;
      case EXPRINFO_NORM:
        if(n == 0)
        {
          s.degree = -1;
          s.property_flags = __s.s.property_flags;
          s.property_flags.c_info = convex_e();
          s.annotation_flags.integer = false;
        }
        else
        { 
          if(s.property_flags.separable != t_true ||
             __s.s.property_flags.separable != t_true)
            s.property_flags.separable = t_maybe;
        }
        if(__s.flags & SIMPLIFY_0_IS_PROD &&
                __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE)
        {
          r.coeffs[n-ndel] *= pow(fabs(__s.m.nd()), r.params.nd());
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.coeffs[n-ndel] *= pow(fabs(__s.m.nd()), r.params.nd());
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        break;
      case EXPRINFO_POW:
        {
          bool spupd = false;
          if(n == 0)
          {
            if(__s.flags & SIMPLIFY_0_IS_CONST)
            { // a^x = exp(x * ln a)
              double help =
                  __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                               __s.m.nd();
              r.operator_type = EXPRINFO_EXP;
              r.coeffs[0] = log(r.coeffs[0] * help) * r.coeffs[1];
              r.coeffs.pop_back();
              (*delnodes).push_back(__s.nnum);
              s.degree = -1;
              ndel++;
            }
            else
              s.degree = __s.s.degree;
            if((__s.flags & SIMPLIFY_0_IS_SUM &&
                 __s.flags & SIMPLIFY_0_SUM_IS_SIMPLE) ||
               (__s.flags & SIMPLIFY_0_IS_PROD &&
                __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE))
            {
              transfer = __s.transfer;
              m = __s.m;
              nnum = __s.nnum;
              flags = __s.flags;
              //simple_sum_prod_update(r, __s);
            }
            else
              transfer = NULL;
          }
          if(n != 0 && r.params.nd() == 0 &&
             __s.flags & SIMPLIFY_0_IS_CONST &&
             __s.flags & SIMPLIFY_0_CONST_IS_INTEGER)
          {
            double help = __s.m.nn()*r.coeffs[n];
            r.coeffs[n] = 1.;
            long int helpl = lrint(help);
            if(helpl - help == 0 && helpl < MAXINT && helpl > -MAXINT)
            {
              switch(helpl)
              {
#if 0
                // There is still a lurking bug
                case 0:
                  r.operator_type = EXPRINFO_CONSTANT;
                  r.params = 1.;
                  r.coeffs.erase(r.coeffs.begin(),r.coeffs.end());
                  (*delnodes).push_back(__s.nnum);
                  s.degree = 0;
                  ndel++;
                  break;
                case 1:
                  (*delnodes).push_back(__s.nnum);
                  (*delnodes).push_back(r.node_num);
                  break;
#endif
                case 2:
                  r.operator_type = EXPRINFO_SQUARE;
                  (*delnodes).push_back(__s.nnum);
                  s.degree = s.degree == -1 ? -1 : s.degree*2;
                  r.coeffs.pop_back();
                  ndel++;
                  spupd = true;
                  break;
                case -1:
                  r.params = r.coeffs[0];
                  r.coeffs[0] = 1.;
                  r.operator_type = EXPRINFO_INVERT;
                  (*delnodes).push_back(__s.nnum);
                  s.degree = s.degree == -1 ? -1 : s.degree*__s.m.nn();
                  flags = SIMPLIFY_0_IS_MULTIPLICATIVE;
                  r.coeffs.pop_back();
                  ndel++;
                  break;
                default:
                  r.params = (int) helpl;
                  r.operator_type = EXPRINFO_INTPOWER;
                  (*delnodes).push_back(__s.nnum);
                  s.degree = s.degree == -1 ? -1 : s.degree*__s.m.nn();
                  r.coeffs.pop_back();
                  ndel++;
                  break;
              }
            }
            //else
            //  spupd = true;
          }
          //else if(n != 0)
          //  spupd = true;
          if(spupd)
            simple_sum_prod_update(r, *this);
          if(n != 0)
          {
            if(flags != SIMPLIFY_0_IS_MULTIPLICATIVE)
              flags = 0;
            nnum = r.node_num;
          }
          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.;
          }
          s.property_flags.c_info = convex_e();
          s.property_flags.separable = t_maybe;
        }
        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);
          s.degree = 0;
        }
        else
          simple_sum_prod_update(r, __s);
        if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]) &&
           is_integer(r.params.nd()))
          s.annotation_flags.integer = true;
        break;
      case EXPRINFO_SQROOT:
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          help *= r.coeffs[0];
          help += r.params.nd();
          r.params = sqrt(help);
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        else
          simple_sum_prod_update(r, __s);
        break;
      case EXPRINFO_ABS:
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          help *= r.coeffs[0];
          help += r.params.nd();
          r.params = fabs(help);
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        else
          simple_sum_prod_update(r, __s);
        if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]) &&
           is_integer(r.params.nd()))
          s.annotation_flags.integer = true;
        break;
      case EXPRINFO_EXP:
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          help *= r.coeffs[0];
          help += r.params.nd();
          r.params = exp(help);
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        else
          simple_sum_prod_update(r, __s);
        break;
      case EXPRINFO_LOG:
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          help *= r.coeffs[0];
          help += r.params.nd();
          r.params = log(help);
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        else
          simple_sum_prod_update(r, __s);
        break;
      case EXPRINFO_SIN:
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          help *= r.coeffs[0];
          help += r.params.nd();
          r.params = sin(help);
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        else
          simple_sum_prod_update(r, __s);
        break;
      case EXPRINFO_COS:
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          help *= r.coeffs[0];
          help += r.params.nd();
          r.params = cos(help);
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        else
          simple_sum_prod_update(r, __s);
        break;
      case EXPRINFO_ATAN2:
        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.;
        }
        // const should be computed here, too.
        break;
      case EXPRINFO_GAUSS:
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          help *= r.coeffs[0];
          help -= r.params.d()[0];
          help /= r.params.d()[1];
          r.params = exp(-help*help);
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        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];
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_PROD &&
                __s.flags & SIMPLIFY_0_PROD_IS_SIMPLE)
        {
          r.coeffs[0] *= __s.m.nd();
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.coeffs[0] *= __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        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);
          s.degree = 0;
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.params = r.params.nd() / __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        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;
          }
        }
        if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          if(n == 0)
            r.params = r.params.nd() * __s.m.nd();
          else
            r.params = r.params.nd() / __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        break;
      case EXPRINFO_INTPOWER:
        if(r.params.nn() > 0)
          s.degree = __s.s.degree == -1 ? -1 : __s.s.degree*r.params.nn();
        else
          s.degree = -1;
        if(__s.flags & SIMPLIFY_0_IS_CONST)
        { // compute constant expressions
          double help =
                __s.flags & SIMPLIFY_0_CONST_IS_INTEGER ? __s.m.nn() :
                                                          __s.m.nd();
          r.operator_type = EXPRINFO_CONSTANT;
          r.params = pow(help, r.params.nn());
          (*delnodes).push_back(__s.nnum);
          s.degree = 0;
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.coeffs[0] *= __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        if(__s.s.annotation_flags.integer && is_integer(r.coeffs[0]))
          s.annotation_flags.integer = true;
        break;
      case EXPRINFO_IN:
        if(n == 0)
          r_idx = -1;
        else
        {
          if(l_nnum == __s.nnum)
          {
            if(r_idx == -1)
              r_idx = n-ndel-1;
            r.params.i()[n-ndel] /= r.coeffs[n-ndel];
            r.coeffs[n-ndel] = 1.;
            r.params.i()[r_idx].intersect(r.params.i()[n-ndel]);
            (*deledges).push_back(std::make_pair(nnum, __s.nnum));
          }
          else if(r_idx != -1)
          {
            int ne_del = n-ndel-r_idx-1;
            r.coeffs.erase(r.coeffs.begin()+r_idx+1,
                           r.coeffs.begin()+n-ndel);
            r.params.i().erase(r.params.i().begin()+r_idx+1,
                               r.params.i().begin()+n-ndel);
            ndel += ne_del;
            r_idx = -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();
          help *= r.coeffs[n-ndel];
          if(r.params.i()[n-ndel].contains(help))
          {
            r.params.i().erase(r.params.i().begin()+n-ndel);
            (*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();
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        else if(__s.flags & SIMPLIFY_0_IS_MULTIPLICATIVE)
        {
          r.coeffs[n-ndel] *= __s.m.nd();
          if(!(__s.flags & SIMPLIFY_0_IS_CORRECTEDMULT))
            __mod->node(__s.nnum)->params = 1.;
        }
        else if(__s.flags & SIMPLIFY_0_IS_SUM &&
           __s.flags & SIMPLIFY_0_SUM_IS_SIMPLE)
        {
          double help = r.coeffs[n-ndel]*__s.m.nd();
          r.params.i()[n-ndel] -= help;
          r.coeffs[n-ndel] *= (*__s.transfer)[0];
          if(__s.flags & SIMPLIFY_0_IS_GHOST)
            transfer_ghost_down(__s.nnum, r.node_num);
          else
            (*delnodes).push_back(__s.nnum);
        }
        if(r.coeffs[n-ndel] != 1)
        {
          r.params.i()[n-ndel] /= r.coeffs[n-ndel];
          r.coeffs[n-ndel] = 1.;
        }
        break;
      case EXPRINFO_POLY:
      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:
        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.;
        }
      case EXPRINFO_EXPECTATION:
      case EXPRINFO_INTEGRAL:
        std::cerr << "simplify_visitor_0: Integral Operations NYI" << std::endl;
        throw "NYI";
        break;
      case EXPRINFO_LOOKUP:
      case EXPRINFO_PWLIN:
      case EXPRINFO_SPLINE:
      case EXPRINFO_PWCONSTLC:
      case EXPRINFO_PWCONSTRC:
        std::cerr << "simplify_visitor_0: Table Operations NYI" << std::endl;
        throw "NYI";
        break;
      case EXPRINFO_DET:
      case EXPRINFO_COND:
      case EXPRINFO_PSD:
      case EXPRINFO_MPROD:
      case EXPRINFO_FEM:
        std::cerr << "simplify_visitor_0: Matrix Operations NYI" << std::endl;
        throw "NYI";
        break;
      case EXPRINFO_RE:
      case EXPRINFO_IM:
      case EXPRINFO_ARG:
      case EXPRINFO_CPLXCONJ:
        std::cerr << "simplify_visitor_0: Complex Numbers NYI" << std::endl;
        throw "NYI";
        break;
      case EXPRINFO_CMPROD:
      case EXPRINFO_CGFEM:
        std::cerr << "simplify_visitor_0: Constant Matrix Operations NYI" << std::endl;
        throw "NYI";
        break;
      case EXPRINFO_LIN:
        break;
      case EXPRINFO_QUAD:
        break;
      default:
        std::cerr << "Simplifier: Unknown Operator Type " << r.operator_type
          << std::endl;
        throw "NYI";
        break;
    }
  }
  else if(r.operator_type > 0)
  {
    ;
  }
  n++;
  v_i.set_all(__s.v_i);
  l_nnum = __s.nnum;
  l_flags = __s.flags;
  l_s = __s.s;
}

bool model::basic_simplify_0()
{
  std::vector<unsigned int> delnodes;
  std::vector<std::pair<unsigned int,unsigned int> > deledges;
  simplify_visitor_0 w(gid->number_of_variables(), &delnodes, &deledges, this);
  bool result = false;

  recursive_postorder_walk(ground(), w);
  result = !delnodes.empty() || !deledges.empty();
  for(std::vector<std::pair<unsigned int,unsigned int> >::iterator b =
      deledges.begin(); b != deledges.end(); ++b)
    remove_edge(node((*b).first),node((*b).second));
  for(std::vector<unsigned int>::iterator b = delnodes.begin(); b != delnodes.end();
      ++b)
    remove_node(*b);
  return result;
}

bool model::basic_simplify_1_merge(walker __s, std::list<walker>& _todo,
                                   std::vector<int>& _combined, std::vector<bool>& _fld,
                                   std::vector<int>& _sorted)
{
  std::list<walker> dew;
  walker _ew, _fw;
  parents_iterator _p, _q;
  bool result = false;
  unsigned int snn = __s == sky() ? number_of_nodes() : __s->node_num;

  if(_sorted[snn] <= 0)
  {
    _sorted[snn] += 2; // 0 -> 2, -1 -> 1
    sort_parent_edges(__s, expression_node::parents_compare());
  }
  for(_p = __s.parent_begin(); _p != __s.parent_end(); ++_p)
  {
    _ew = __s << _p;
    if(_combined[_ew->node_num] != (int)_ew->n_children)
      continue;
    if(abs(_sorted[_ew->node_num]) != 1)
    {
      _sorted[_ew->node_num]--; // 2 -> 1, 0 -> -1
      if(_ew->operator_type == EXPRINFO_SUM ||
         _ew->operator_type == EXPRINFO_PROD ||
         _ew->operator_type == EXPRINFO_MAX ||
         _ew->operator_type == EXPRINFO_MIN ||
         _ew->operator_type == EXPRINFO_MONOME ||
         _ew->operator_type == EXPRINFO_NORM ||
         _ew->operator_type == EXPRINFO_MEAN)
      {
        // this is an awful hack, but I have not found a better
        // one. The cleanest way would be to fully work out the
        // "real" dag class which depends on an edge class
        std::map<std::pair<unsigned int, unsigned int>,
                 std::pair<double, int> > __nr_c_ref;
        std::map<std::pair<unsigned int, unsigned int>,
                 std::pair<double, int> >::iterator __f;
        children_iterator _c;
        int n;
        
        for(_c = _ew.child_begin(), n=0; _c != _ew.child_end(); ++_c, ++n)
        {
          unsigned int cnn = (_ew >> _c)->node_num, _h;
          for(_h = 0; ; _h++)
          {
            __f = __nr_c_ref.find(std::make_pair(cnn, _h));
            if(__f == __nr_c_ref.end())
              break;
          }
          if(_ew->operator_type == EXPRINFO_MONOME)
            __nr_c_ref.insert(std::make_pair(std::make_pair(cnn, _h),
                           std::make_pair(_ew->coeffs[n],_ew->params.n()[n])));
          else
            __nr_c_ref.insert(std::make_pair(std::make_pair(cnn, _h),
                              std::make_pair(_ew->coeffs[n],0)));
        }
        sort_child_edges(_ew, expression_node::children_compare());
        for(_c = _ew.child_begin(), n=0; _c != _ew.child_end(); ++_c, ++n)
        {
          unsigned int cnn = (_ew >> _c)->node_num, _h;

          for(_h = 0; ; _h++)
          {
            __f = __nr_c_ref.find(std::make_pair(cnn, _h));
#if 0
          if(__f == __nr_c_ref.end())
          {
            std::cerr << "FATAL Error: This can't happen ;) in basic_simplify_1_merge"
              << std::endl;
            exit(-42);    // guess why!
          }
#endif
            if(__f != __nr_c_ref.end())
              break;
          }
          _ew->coeffs[n] = __f->second.first;
          if(_ew->operator_type == EXPRINFO_MONOME)
            _ew->params.n()[n] = __f->second.second;
          __nr_c_ref.erase(__f);
        }
      }
    }
  }
  for(_p = __s.parent_begin(); _p != __s.parent_end(); ++_p)
  {
    _ew = __s << _p;
    if(_combined[_ew->node_num] == (int)_ew->n_children)
      break;
  }
  while(_p != __s.parent_end())
  {
    for(_q = __s.parent_begin(); _q != __s.parent_end(); ++_q)
    {
      _ew = __s << _q;
      if(_combined[_ew->node_num] == (int)_ew->n_children)
        break;
    }
    for(++_p; _p != __s.parent_end(); ++_p)
    {
      _fw = __s << _p;
      if(_combined[_fw->node_num] != (int)_fw->n_children)
        continue;
      if(_ew.node() == _fw.node())
      {
        _q = _p;
        _ew = __s << _q;
        continue;
      }
      if(expression_node::parents_compare_eq()(*_ew, *_fw) &&
         children_eq(_ew, _fw))
        // up until now they are all equal
        dew.insert(dew.end(), _fw);
      else
      { // now they have changed
        if(dew.size() != 0)
          /* more than one have to be merged */
          break;
        else
        {
          /* no merging necessary, so flood it unless it's already flooded */
          if(!_fld[_ew->node_num])
          {
            _fld[_ew->node_num] = true;
            basic_simplify_1_flood(_ew, _combined);
            if(_ew->operator_type != EXPRINFO_CONSTANT)
              _todo.push_back(_ew);
          }
          _q = _p;
          _ew = __s << _q;
        }
      }
    }
    if(dew.size() == 0)
    { /* we are at the end of the parents list */
      if(!_fld[_ew->node_num])
      {
        _fld[_ew->node_num] = true;
        basic_simplify_1_flood(_ew, _combined);
        if(_ew->operator_type != EXPRINFO_CONSTANT)
          _todo.push_back(_ew);
      }
      break;
    }

    /* now merge */
    _ew->n_parents += dew.size();
    for(std::list<walker>::iterator _b = dew.begin(); _b != dew.end(); ++_b)
    {
      std::list<walker>::iterator _hlp;
      
      free_node_num((*_b)->node_num);

      for(_hlp = _todo.begin(); _hlp != _todo.end(); ++_hlp)
      {
        if(*_hlp == *_b)
        {
          _todo.erase(_hlp);
          break;
        }
      }

      if(objective == *_b)
        objective = _ew;
      
      // make sure that constraint numbering still works
      unsigned int cn;
      if(get_const_num((*_b)->node_num, cn))
      {
        unsigned int cc;
        if(get_const_num(_ew->node_num, cc))
        {
          gid->unused_constr(const_name(cn));
          gid->remove_const_id(cn);
        }
        else
          gid->make_const_back_ref(_ew->node_num, cn);
        // remove the constraint from the constraints array
        for(std::vector<walker>::iterator __x = constraints.begin();
            __x != constraints.end(); ++__x)
          if(*__x == *_b)
          {
            constraints.erase(__x);
            break;
          }
      }
      
#if MODEL_INLINE_DEBUG_SIMPLIFY
      std::cout << "Merging " << _ew->node_num << " and " <<
        (*_b)->node_num << std::endl;
#endif
      merge(_ew, *_b, false);
      result = true;
    }
    dew.erase(dew.begin(), dew.end());
    if(_ew->operator_type == EXPRINFO_VARIABLE)
      store_variable(_ew);
    else if(_ew->operator_type == EXPRINFO_GHOST)
      store_ghost(_ew);

    /* this one is finished, so flood it unless it's already flooded */
    if(!_fld[_ew->node_num])
    {
      _fld[_ew->node_num] = true;
      basic_simplify_1_flood(_ew, _combined);
      if(_ew->operator_type != EXPRINFO_CONSTANT)
        _todo.push_back(_ew);
    }
    for(_p = __s.parent_begin(); _p != __s.parent_end(); ++_p)
    {
      _fw = __s << _p;
      if(_combined[_fw->node_num] == (int)_fw->n_children)
        break;
    }
  }
  return result;
}

bool model::basic_simplify_1()
{
  std::vector<int> _combined(number_of_nodes(), 0);
  std::vector<bool> _fld(number_of_nodes(), false);
  std::vector<int> _sorted(number_of_nodes()+1, 0);
  std::list<walker> _todo(1, sky());
  std::list<walker>::iterator __xxcl, __help;
  bool result = false;

  while(!_todo.empty())
  {
    for(__xxcl = _todo.begin(); __xxcl != _todo.end(); __xxcl = __help)
    {
#if MODEL_INLINE_DEBUG_SIMPLIFY
      for(std::list<walker>::const_iterator __rr = _todo.begin();
          __rr != _todo.end(); ++__rr)
        std::cerr << (*__rr)->node_num << ", ";
      std::cerr << std::endl;
#endif
      __help = __xxcl;
      ++__help;
      if((*__xxcl).is_ground()||(*__xxcl).is_root())
        // everything ends at root nodes or above
        _todo.erase(__xxcl);
      else if((*__xxcl).is_sky() ||
              basic_simplify_1_is_ready(*__xxcl, _combined, false))
      { /* some of the parents of this node are ready to be merged */
        result = basic_simplify_1_merge(*__xxcl, _todo, _combined, _fld,
                                        _sorted) || result;
        if(basic_simplify_1_is_ready(*__xxcl, _combined, true))
        /* all of the parents of this node are merged */
          _todo.erase(__xxcl);
      }
    }
  }
  return result;
}

bool model::basic_simplify_2()
{
  //walker h;

  return false;
}

void model::prepare_after_read()
{
  walker _w;
  unsigned int cnum;

  compress_numbers();
  for(ref_iterator _x = node_ref.begin(); _x != node_ref.end(); ++_x)
  {
    _w = *_x;
    if(_w->f_bounds.inf() != -INFINITY || _w->f_bounds.sup() != INFINITY)
    {
      constraints.insert(constraints.end(), _w);
      if(!gid->get_const_num(_w->node_num, cnum))
        cnum = next_constraint_num();
      gid->mk_gconstref(cnum, _w);
    }
  }
  
  arrange_constraints();
}

bool model::basic_simplify()
{
  bool have_simplified = false, first = true;
  int have_simp;

  have_simp = 7;
  do
  {
    have_simp &= ~1;
    if((first || have_simp) && basic_simplify_0())
      have_simp |= 1;
#if MODEL_INLINE_DEBUG_SIMPLIFY
    write();
#endif
    have_simp &= ~2;
    if((first || have_simp) && basic_simplify_1())
      have_simp |= 2;
#if MODEL_INLINE_DEBUG_SIMPLIFY
    write();
#endif
    have_simp &= ~4;
    if((first || have_simp) && basic_simplify_2())
      have_simp |= 4;
#if MODEL_INLINE_DEBUG_SIMPLIFY && 0
    write();
#endif
    have_simplified = have_simplified || (have_simp != 0);
    first = false;
  } while(have_simp != 0);

  return have_simplified;
}

void model::write(std::ostream &o) const
{
  std::vector<bool> _nprinted(number_of_nodes(), false);
  expression_print_visitor p(_nprinted, o);
  int precsave = o.precision(18);
  
  o << "<N> V " << gid->number_of_variables() << std::endl;
  for(children_iterator b = ground().child_begin(); b != ground().child_end();
      b++)
  {
    recursive_preorder_walk_if(ground() >> b, p);
  }
  o << "<N> M ";
  if((const_walker) objective != ground())
    o << objective->node_num;
  else
    o << '0';
  o << ' ';
  if(ocoeff == -1)
    o << "max";
  else if(ocoeff == 1)
    o << "min";
  else
    o << "m00";
  o << std::endl;

  for(unsigned int i=0; i < number_of_variables(); ++i)
    o << "<N> N " << i << " '" << gid->var_name(i) << "'" << std::endl;
  for(unsigned int i=0; i < gid->n_fixed_vars(); ++i)
    o << "<N> F '" << gid->fixed_var(i).first << "' "
      << gid->fixed_var(i).second << std::endl;
  for(unsigned int i=0; i < gid->n_unused_vars(); ++i)
    o << "<N> U '" << gid->unused_var(i) << "'" << std::endl;
  for(unsigned int i=0; i < gid->n_unused_constrs(); ++i)
    o << "<N> S '" << gid->unused_constr(i) << "'" << std::endl;
  for(std::vector<walker>::const_iterator __b = constraints.begin();
      __b != constraints.end(); ++__b)
  {
    unsigned int cn;
    if(!gid->get_const_num((*__b)->node_num, cn))
    {
      std::cerr << "FATAL in write: Could not get constraint number for node "
        << (*__b)->node_num << std::endl;
      throw "Programming error";
    }
    o << "<N> c " << cn << " " << (*__b)->node_num << std::endl;
    o << "<N> C " << cn << " '" << const_name(cn) << "'" << std::endl;
  }
#if 0
  for(unsigned int i=0; i < gid->number_of_constraints(); ++i)
  {
    if(constraint(i) == (const walker&) ground()) continue;
    o << "<N> c " << i << " " << constraint(i)->node_num << std::endl;
    o << "<N> C " << i << " '" << const_name(i) << "'" << std::endl;
  }
#endif
  if(ocoeff != 0 || gid->obj_adj() != 0)
  {
    if(ocoeff == 0) gid->obj_mult(1);
    o << "<N> O '" << gid->obj_name() << "'";
    if(gid->obj_adj() != 0 || gid->obj_mult() != 1)
      o << ": d " << gid->obj_adj();
    if(gid->obj_mult() != 1)
      o << ", " << gid->obj_mult();
    o << std::endl;
  }
  o.precision(precsave);
}

model::walker model::variable(unsigned int _vnum)
{
  bool make_new(false);
  walker __r(ground());
  ref_iterator __x(var_ref.end());

  if(!gid->have_gvar_ref(_vnum))
  {
    new_variables(_vnum+1);
    make_new = true;
  }
  else
  {
    __x = lower_bound(var_ref.begin(), var_ref.end(), _vnum,
                      __docompare_variables());
    if(__x != var_ref.end() && (*__x)->params.nn() == (int)_vnum)
    {
      __r = *__x;
      make_new = false;
    }
    else if(gid->empty(var(_vnum)))
      make_new = true;
  }
  if(make_new)
  {
    expression_node _en(EXPRINFO_VARIABLE, next_num());
    _en.v_i.reserve(gid->number_of_variables());
    _en.v_i.set(_vnum);
    _en.params = (int)_vnum;
    __r = between_back(ground(), sky(), _en);
    var_ref.insert(__x, __r);
    return store_node(__r);
  }
  else if(__r == ground()) // we have to produce a ghost variable
  {
    __r = var(_vnum);
    __r = ghost(__r->node_num);
    // var_ref.insert(__x, __r); removed Oct. 20 (ghosts must not be in var_ref)
    return __r;
  }
  else
    return __r;
}

model::walker model::vnary(int expr_type, ...)
{
  va_list _a_list;
  expression_node _en(expr_type, next_num());
  std::vector<walker> _args;
  const walker* _e;
  int __i;

  va_start(_a_list, expr_type);
  for(__i=0;;__i++)
  {
    _e = va_arg(_a_list, const walker*);
    if(_e == EXPR_LASTARG) break;
    _args.push_back(*_e);
    _en.v_i.set_all((*_e)->v_i);
  }
  va_end(_a_list);
  _en.coeffs.insert(_en.coeffs.end(),__i,1.0);
  if(expr_type == EXPRINFO_PROD)
    _en.params = 1.;
  else
    _en.params = 0.;
  return store_node(split_back(ground(), _args, _en));
}

void model::detect_0chain()
{
  unsigned int f = 1;
  unsigned int n = number_of_nodes();
  std::vector<unsigned int> t(n,0);
  std::vector<unsigned int> depth(n,0);
  std::vector<unsigned int> base;

  detect_0chain_visitor d(t, depth, base, f);

  evaluate(d, ground());
  for(unsigned int i = 0; i < n; ++i)
  {
    if(gid->empty(gid->node(i)))
      continue;
    if(t[i] == 0)
    {
      gid->node(i)->sem.info_flags.has_0chnbase = t_false;
      gid->node(i)->sem._0chnbase = i;
      continue;
    }
    if(base[t[i]-1] == i)
      gid->node(i)->sem.info_flags.has_0chnbase = t_false;
    else
      gid->node(i)->sem.info_flags.has_0chnbase = 
                                        (depth[i] == 1 ? t_maybe : t_true);
    gid->node(i)->sem._0chnbase = base[t[i]-1];
  }
}

---