standardize.cpp

Go to the documentation of this file.

/* $Id: standardize.cpp 1147 2015-05-04 14:01:51Z stefan $
 *
 * Name:    standardize.cpp
 * Author:  Pietro Belotti
 * Purpose: standardize all expressions in a problem
 *
 * (C) Carnegie-Mellon University, 2006-11.
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include <vector>

#include "CoinHelperFunctions.hpp"

#include "BonBabSetupBase.hpp"

#include "CouenneTypes.hpp"
#include "CouenneExpression.hpp"
#include "CouenneExprIVar.hpp"
#include "CouenneExprSub.hpp"
#include "CouenneExprClone.hpp"
#include "CouenneProblem.hpp"
#include "CouenneProblemElem.hpp"
#include "CouenneDepGraph.hpp"
#include "CouenneSdpCuts.hpp"

using namespace Ipopt;
using namespace Couenne;

void replace (CouenneProblem *p, int wind, int xind);


bool CouenneProblem::standardize () {

  if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
    printf ("Reformulation. current point: %d vars -------------------\n", 
            domain_.current () -> Dimension ());
    for (int i=0; i<domain_.current () -> Dimension (); i++)
      printf ("%3d %20g [%20g %20g]\n", i, domain_.x (i), domain_.lb (i), domain_.ub (i));
  }

  bool retval = true;

  // create dependence graph to assign an order to the evaluation (and
  // bound propagation, and -- in reverse direction -- implied bounds)
  graph_ = new DepGraph;

  for (std::vector <exprVar *>::iterator i = variables_ . begin ();
       i != variables_ . end (); ++i)
    graph_ -> insert (*i);

  // allocate space in auxiliaries_ from commonexprs_

  //int initVar = variables_ . size () - commonexprs_ . size ();

  //std::set <int> DefVarDiffSet;

  // // DEFINED VARIABLES -----------------------------------------------------------------------

  // // standardize initial aux variables (aka defined variables, aka
  // // common expression)

  // if (commonexprs_.size ()) jnlst_ -> Printf (J_ALL, J_REFORMULATE,
  //                                          "%d common exprs, initVar = %d = %d - %d\n", 
  //                                          commonexprs_ . size (), 
  //                                          initVar, 
  //                                          variables_   . size (), 
  //                                          commonexprs_ . size ());

  // for (std::vector <expression *>::iterator i = commonexprs_ . begin ();
  //      i != commonexprs_ . end (); ++i) {

  //   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
  //     printf ("\n=====> [1] stdz common expr [%d] := ", initVar); fflush (stdout);
  //     (*i) -> print (); printf ("\n");
  //   }

  //   exprAux    *naux = (*i) -> standardize (this, false);
  //   expression *img  = (*i);

  //   if (naux)
  //     img = naux -> Image ();

  //   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
  //     printf ("\n=====> [2] "); fflush (stdout);
  //     if (*i)   (*i) -> print (); else printf ("(null)"); printf (" [*i] ");   fflush (stdout);
  //     if (naux) naux -> print (); else printf ("(null)"); printf (" [naux] "); fflush (stdout);
  //     if (img)  img  -> print (); else printf ("(null)"); printf (" [img]\n");
  //   }

  //   // trick to obtain same index as common expression: create exprAux
  //   // with index initVar and replace occurrences with address of
  //   // newly created exprAux through auxiliarize()

  //   exprVar *newvar;
  //   //img -> isInteger () ? exprAux::Integer : exprAux::Continuous);

  //   //auxiliarize (newvar); // puts newvar at right position in variables_

  //   if (((*i) -> Type () == VAR) ||
  //    ((*i) -> Type () == AUX)) {

  //     newvar = new exprAux (img, initVar, 1 + img -> rank (), exprAux::Unset, &domain_);
  //     replace (this, initVar, img -> Index ());
  //     //auxiliarize (variables_ [initVar], 
  //     //variables_ [img -> Index ()]);

  //     //delete variables_ [initVar];

  //     variables_ [initVar] = newvar;
  //     variables_ [initVar] -> zeroMult ();

  //   } else {

  //     if (img -> dependsOn (&initVar, 1, TAG_AND_RECURSIVE)) {

  //    //printf ("depends! "); img -> print (); 

  //    expression *diff = new exprSub (new exprClone (variables_ [initVar]), img);

  //    //printf ("; diff = "); diff -> print ();

  //    exprAux *diffAux = diff -> standardize (this, false);

  //    //printf ("; aux: "); if (diffAux) diffAux -> print (); 

  //    //if (diffAux)
  //    exprAux *newAux = addAuxiliary (diff);

  //    //printf ("; real aux: "); if (newAux) newAux -> print (); putchar ('\n');
        
  //    //Lb (newAux -> Index ()) = 
  //    //Ub (newAux -> Index ()) = 0.;

  //    DefVarDiffSet. insert (newAux -> Index ());

  //     } else {

  //    newvar = new exprAux (img, initVar, 1 + img -> rank (), exprAux::Unset, &domain_);
  //    //replace (this, initVar, newvar -> Index ());

  //    auxiliarize (newvar);

  //    //delete variables_ [initVar];
  //    variables_ [initVar] = newvar;

  //    graph_ -> insert (newvar);
  //    //if (naux) 
  //    graph_ -> erase (naux);
  //     }
  //   }

  //   //variables_ . erase (variables_ . end () - 1);

  //   if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
  //     if (naux) {
  //    printf ("done: "); fflush (stdout);
  //    naux -> print ();
  //    printf (" := "); fflush (stdout);
  //    naux -> Image () -> print (); printf ("\n..."); fflush (stdout);
  //     } else if (*i) {
  //    (*i) -> print ();
  //    //printf (" := "); fflush (stdout);
  //    //aux -> Image () -> print (); 
  //    printf ("\n");
  //     } else printf ("[n]aux NULL!\n");
  //   }

  //   initVar++;
  // }

  // OBJECTIVES ------------------------------------------------------------------------------

  for (std::vector <CouenneObjective *>::iterator i = objectives_.begin ();
       i != objectives_.end (); ++i) {

    if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
      printf ("Objective [code: %d]", (*i) -> Body () -> code ()); 
      (*i) -> print ();
    }

    // check for constant objective

    std::set <int> deplist;

    if (0 == (*i) -> Body () -> DepList (deplist,  TAG_AND_RECURSIVE)) {

      // This objective is constant. Store the value and use it in
      // CouenneSolverInterface.

      constObjVal_ = (*i) -> Body () -> Value ();

    } else {

      exprAux *aux = (*i) -> standardize (this);

      if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
        printf (" objective "); (*i) -> print (); 
        if (aux) {printf (" admits aux "); aux -> print ();}
      }

      if (aux) {
        //delete ((*i) -> Body ());
        (*i) -> Body (new exprClone (aux));
      }

      if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
        printf (". New obj: "); (*i) -> print (); printf ("\n");
      }
    }
  }

  // CONSTRAINTS ----------------------------------------------------------------------------

  // commuted_ is an array with a flag for each original variable,
  // which is true at position i if initially original variable x_i
  // became auxiliary (because of constraint 
  // 
  // x_{k+1} + f(x_1,x_2...,x_k} <=/=/>= 0
  //
  // becoming
  //
  // x_{k+1} <=/=/>= - f(x_1,x_2...,x_k}

  commuted_ = new bool [nVars ()];
  CoinFillN (commuted_, nVars (), false);

  std::vector <std::vector <CouenneConstraint *>::iterator> iters2erase;

  for (std::vector <CouenneConstraint *>::iterator i = constraints_.begin (); 
       i != constraints_.end (); ++i) {

    if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
      printf ("\nReformulating constraint: "); 
      (*i) -> print ();
    }

    CouNumber 
      conLb = (*((*i) -> Lb ())) (),
      conUb = (*((*i) -> Ub ())) ();

    // sanity check: if (due to bad model or redundancies) the
    // constraint's body is constant, delete it -- check if it's
    // within bounds.

    expression *eBody = (*i) -> Body ();

    if (eBody -> Linearity () <= CONSTANT) {

      CouNumber bodyVal = (*eBody)();

      if ((bodyVal < conLb - COUENNE_BOUND_PREC) ||
          (bodyVal > conUb + COUENNE_BOUND_PREC)) { // all variables eliminated, but out of bounds
        
        jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, 
                          "Constraint: all variables eliminated, but value %g out of bounds [%g,%g]: ",                           
                          bodyVal, conLb, conUb);

        if (jnlst_ -> ProduceOutput (J_SUMMARY, J_PROBLEM))
          (*i) -> print ();

        retval = false;
        break;

      } else {

        iters2erase.push_back (i);
        continue; // all variables eliminated and constraint is redundant
      }
    }

    exprAux *aux = (*i) -> standardize (this);

    if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
      printf (" reformulated: aux w[%d] ", aux ? (aux -> Index ()) : -2); 
      (*i) -> print ();
    }

    if (aux) { // save if standardized

      // this is a top level auxiliary, i.e., an auxiliary whose node
      // in the DAG stands at the maximum level -- no other variable
      // depends on it as it is the lhs of a constraint.

      aux -> top_level () = true;

      // this constraint f(x)<=b is therefore replaced by w<=b, where
      // w is the top-level auxiliary associated with f(x)

      //printf ("delete %x: ", ((*i) -> Body ())); ((*i) -> Body ()) -> print ();
      //printf ("\n"); 
      //delete ((*i) -> Body ());

      (*i) -> Body (new exprClone (aux));

      //      con2.push_back (*i);
    }
    else {

      // left-hand side not reformulated, therefore this is probably
      // an affine expression

      //CouNumber lb, ub;

      //printf ("let's see what's left of the body: "); fflush (stdout);
      //(*i) -> Body () -> print (); printf ("\n");

      // (*i) -> Body () -> getBounds (lb, ub);

      // if ((((*((*i) -> Lb ())) ()) > ub + COUENNE_EPS) ||
      //          (((*((*i) -> Ub ())) ()) < lb - COUENNE_EPS)) {

      //        jnlst_ -> Printf (J_SUMMARY, J_PROBLEM, "found infeasible constraint [%g,%g] vs [%g,%g]\n", 
      //                          lb, ub,
      //                          ((*((*i) -> Lb ())) ()),
      //                          ((*((*i) -> Ub ())) ()));

      //        if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE))
      //          (*i) -> print ();

      //        retval = false;
      // }

      // delete constraint as no aux associated with it. It is a
      // linear constraint now (after reformulation), so it will be
      // inserted as a linearization of the aux's definition (an
      // affine function).

      iters2erase.push_back (i);
    }

    //(*i) -> Body () -> realign (this);

    if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {printf (" --> "); (*i) -> print (); printf ("\n\n");}

    /*printf ("=== "); fflush (stdout); 
    aux -> print (); printf (" := "); fflush (stdout);
    aux -> Image () -> print (); printf ("\n");*/
  }

  for (unsigned int i = iters2erase.size (); i--;)
    constraints_. erase (iters2erase [i]);

  if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) { 
    // Use with caution. Bounds on auxs are not defined yet, so valgrind complains
    printf ("Done with standardization (careful, bounds below can be nonsense):\n");
    print (); 
  }

  // Extra auxiliaries: useless without CouenneSDPcuts

  if (bonBase_) {

    int freq;

    if ((bonBase_ -> options () -> GetIntegerValue ("sdp_cuts", freq, "couenne.")) &&
        (freq != 0))

      sdpCutGen_ = new CouenneSdpCuts (this, jnlst_, bonBase_ -> options ());
  }

  // TODO: move redundancy elimination here. Then add another round of simplification

#define ANTICIPATE_DELETE_REDUND
#ifdef  ANTICIPATE_DELETE_REDUND

  bool has_changed;

  do {

    has_changed = false;

    // TODO: resolve duplicate index in exprQuad before restoring this

    // remove duplicates

    std::string delete_redund;

    if (bonBase_) bonBase_ -> options () -> GetStringValue ("delete_redundant", delete_redund, "couenne."); 
    else  delete_redund = "yes";

    if (delete_redund == "yes")

      // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
      for (std::vector <exprVar *>::iterator i = variables_.begin (); 
           i != variables_.end (); ++i)

        if (((*i) -> Multiplicity () > 0) && ((*i) -> Type () == AUX) && ((*i) -> sign () == expression::AUX_EQ)) {

          int type = (*i) -> Image () -> Type ();

          if ((type == VAR) || (type == AUX)) {

            // found w_k = x_h. 
            // 
            // Check if either is integer, the survivor will be integer too
            // Replace all occurrences of w_k with x_h

            /*printf ("redundancy: "); 
              (*i)             -> print (); printf (" := "); 
              (*i) -> Image () -> print (); printf ("\n");*/

            // use the info on the variable to be discarded: tighter
            // bounds and integrality that the replacement might not have.

            int 
              indStays  = (*i) -> Image () -> Index (), // index h
              indLeaves = (*i)             -> Index (); // index k

            if (indStays == indLeaves)  // very strange case, w_h = x_h
              continue;

            // do not swap! x_h could be in turn an auxiliary...
            //
            //if (indStays > indLeaves) 
            //{int swap = indStays; indStays = indLeaves; indLeaves = swap;} // swap

            exprVar 
              *varStays  = variables_ [indStays],
              *varLeaves = variables_ [indLeaves];

            // intersect features of the two variables (integrality, bounds)

            varStays -> lb () = varLeaves -> lb () = CoinMax (varStays -> lb (), varLeaves -> lb ());
            varStays -> ub () = varLeaves -> ub () = CoinMin (varStays -> ub (), varLeaves -> ub ());

            if (varStays  -> isInteger () ||
                varLeaves -> isInteger ()) {

              varStays -> lb () = ceil  (varStays -> lb ());
              varStays -> ub () = floor (varStays -> ub ());

              if (varStays -> Type () == AUX)
                varStays -> setInteger (true);
              else {
                //expression *old = varStays; // !!! leak
                variables_ [indStays] = varStays = new exprIVar (indStays, &domain_);
                auxiliarize (varStays); // replace it everywhere in the problem
                //delete old;
              }
            }

            auxiliarize (varLeaves, varStays); // now replace occurrences of w_k with x_h

            //if (varLeaves -> Index () >= nOrigVars_) // why check? It's not there anymore.
            varLeaves -> zeroMult (); // disable this variable
          }
        }

    // realignment and cleanup of variable might have introduced
    // variable duplication, e.g. x0*x0. Do one round of
    // simplification and reformulate what is simplified.

    for (std::vector <exprVar *>::iterator i = variables_. begin ();
         i != variables_. end (); ++i) 

      if (((*i) -> Multiplicity () > 0) &&
          ((*i) -> Type () == AUX)) {

        expression
          *image = (*i) -> Image (),
          *simpl = image -> simplify ();

        exprVar *subst = NULL;

        if (simpl &&
            (subst = simpl -> standardize (this))) {
          if ((subst -> Type () == VAR) ||
              (subst -> Type () == AUX))
            simpl = new exprClone (subst);
          else simpl = subst;
        }

        if (simpl && simpl != subst) {

          has_changed = true;

          delete (*i) -> Image ();
          (*i) -> Image (simpl);
        }
      }

  } while (has_changed); // should repeat this while there have been
                         // simplifications or duplicates have been
                         // cut
#endif

  //  delete auxSet_;

  // Create evaluation order ////////////////////////////////////////////////////

  // reallocate space for enlarged set of variables
  domain_.current () -> resize (nVars ());

  //graph_ -> print ();
  graph_ -> createOrder ();
  //graph_ -> print ();

  assert (graph_ -> checkCycles () == false);

  // Fill numbering structure /////////////////////////////////////////////////

  int n = nVars ();
  numbering_ = new int [n];
  std::set <DepNode *, compNode> vertices = graph_ -> Vertices ();

  for (std::set <DepNode *, compNode>::iterator i = vertices.begin ();
       i != vertices.end (); ++i)

    numbering_ [(*i) -> Order ()] = (*i) -> Index (); 


  // do initial bound propagation

  for (int i = 0; i < nVars (); i++) {

    int ord = numbering_ [i];

    if (variables_ [ord] -> Multiplicity () <= 0)
      continue;

    if (variables_ [ord] -> Type () == AUX) {

      // initial auxiliary bounds are infinite (they are later changed
      // through branching)

      if (variables_ [ord] -> Index () >= nOrigVars_) { // and one that was not an original, originally...

        domain_.lb (ord) = -COIN_DBL_MAX;
        domain_.ub (ord) =  COIN_DBL_MAX;
      }

      //printf ("from "); variables_ [ord] -> Lb    () -> print (); 

      // tighten them with propagated bounds
      variables_ [ord] -> crossBounds ();

      //printf ("to "); variables_ [ord] -> Lb    () -> print (); printf (", now eval\n");

      if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
        printf (":::: %3d %10g [%10g, %10g]", 
                ord, domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
      }

      // and evaluate them
      domain_.x  (ord) = (*(variables_ [ord] -> Image ())) ();
      domain_.lb (ord) = (*(variables_ [ord] -> Lb    ())) ();
      domain_.ub (ord) = (*(variables_ [ord] -> Ub    ())) ();

      // if (DefVarDiffSet.find (ord) != DefVarDiffSet.end ()) {

      //        domain_.lb (ord) =
      //        domain_.ub (ord) = 0.;
      // }

      if (jnlst_ -> ProduceOutput (J_ALL, J_REFORMULATE)) {
        printf (" --> %10g [%10g, %10g] [", 
                domain_.x (ord), domain_.lb (ord), domain_.ub (ord));
        variables_ [ord] -> Lb () -> print (); printf (",");
        variables_ [ord] -> Ub () -> print (); printf ("]\n");
      }

      bool integer = variables_ [ord] -> isInteger ();

      if (integer) {
        domain_.lb (ord) = ceil  (domain_.lb (ord) - COUENNE_EPS);
        domain_.ub (ord) = floor (domain_.ub (ord) + COUENNE_EPS);
      }
    }
  }

#ifndef ANTICIPATE_DELETE_REDUND

  // TODO: resolve duplicate index in exprQuad before restoring this

  // remove duplicates

  std::string delete_redund;

  if (bonBase_) bonBase_ -> options () -> GetStringValue ("delete_redundant", delete_redund, "couenne."); 
  else  delete_redund = "yes";

  if (delete_redund == "yes")

    // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
    for (std::vector <exprVar *>::iterator i = variables_.begin (); 
       i != variables_.end (); ++i)

    if (((*i) -> Type () == AUX) && ((*i) -> sign () == expression::AUX_EQ)) {

      int type = (*i) -> Image () -> Type ();

      if ((type == VAR) || (type == AUX)) {

        // found w_k = x_h. 
        // 
        // Check if either is integer, the survivor will be integer too
        // Replace all occurrences of w_k with x_h

        /*printf ("redundancy: "); 
        (*i)             -> print (); printf (" := "); 
        (*i) -> Image () -> print (); printf ("\n");*/

        // use the info on the variable to be discarded: tighter
        // bounds and integrality that the replacement might not have.

        int 
          indStays  = (*i) -> Image () -> Index (), // index h
          indLeaves = (*i)             -> Index (); // index k

        if (indStays == indLeaves)  // very strange case, w_h = x_h
          continue;

        // do not swap! x_h could be in turn an auxiliary...
        //
        //if (indStays > indLeaves) 
        //{int swap = indStays; indStays = indLeaves; indLeaves = swap;} // swap

        exprVar 
          *varStays  = variables_ [indStays],
          *varLeaves = variables_ [indLeaves];

        // intersect features of the two variables (integrality, bounds)

        varStays -> lb () = varLeaves -> lb () = CoinMax (varStays -> lb (), varLeaves -> lb ());
        varStays -> ub () = varLeaves -> ub () = CoinMin (varStays -> ub (), varLeaves -> ub ());

        if (varStays  -> isInteger () ||
            varLeaves -> isInteger ()) {

          varStays -> lb () = ceil  (varStays -> lb ());
          varStays -> ub () = floor (varStays -> ub ());

          if (varStays -> Type () == AUX)
            varStays -> setInteger (true);
          else {
            //expression *old = varStays; // !!! leak
            variables_ [indStays] = varStays = new exprIVar (indStays, &domain_);
            auxiliarize (varStays); // replace it everywhere in the problem
            //delete old;
          }
        }

        auxiliarize (varLeaves, varStays); // now replace occurrences of w_k with x_h

        //if (varLeaves -> Index () >= nOrigVars_) // why check? It's not there anymore.
        varLeaves -> zeroMult (); // disable this variable
      }
    }

#endif

  for (int ii=0; ii < nVars (); ii++) {

    int i = numbering_ [ii];

    if ((Var (i) -> Multiplicity () > 0) &&
        (Var (i) -> Type () == AUX) &&
        (Var (i) -> Image () -> isInteger ()) &&
        (Var (i) -> sign () == expression::AUX_EQ))
      Var (i) -> setInteger (true);

    //if (Var (i) -> Multiplicity () == 0)
    //Lb (i) = Ub (i) = 0.;
  }

  // check how many multiplications there are 

//   int nmul = 0;
//   // Look for auxiliaries of the form w:=x and replace each occurrence of w with x
//   for (std::vector <exprVar *>::iterator i = variables_.begin (); 
//        i != variables_.end (); ++i) {

//     if ((*i) -> Type () != AUX ||
//      (*i) -> Multiplicity () <= 0) 
//       continue;

//     expression *img = (*i) -> Image ();

//     if (img -> code () != COU_EXPRMUL) 
//       continue;

//     expression **args = img -> ArgList ();

//     if ((args [0] -> Type () == AUX ||
//       args [0] -> Type () == VAR) &&
//      (args [1] -> Type () == AUX ||
//       args [1] -> Type () == VAR))
//       nmul++;
//   }

//   printf ("\nMULS: %d/%d\n", nmul, variables_.size());

//   exit (-1);

  // TODO: re-compute ranks

  delete [] commuted_;  commuted_ = NULL;
  delete    graph_;     graph_    = NULL;

  return retval;
}