drawCuts.cpp

Go to the documentation of this file.

/* $Id: drawCuts.cpp 490 2011-01-14 16:07:12Z pbelotti $
 *
 * Name:    drawCuts.cpp
 * Author:  Pietro Belotti
 * Purpose: print function and cut information to be later displayed
 *          through PlotUtils' graph
 *
 * (C) Carnegie-Mellon University, 2006. 
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "CouenneCutGenerator.hpp"
#include "CouenneProblem.hpp"
#include "CouenneExprAux.hpp"

using namespace Couenne;

// [cool!] print graph-readable output for displaying inequalities
// on a Cartesian plane

void draw_cuts (OsiCuts &cs, const CouenneCutGenerator *cg, int j, expression *w, expression *img) {

  static bool first_draw = true;
  static CouNumber maxY = -COUENNE_INFINITY,
                   minY =  COUENNE_INFINITY;

  if ((img -> code () == COU_EXPRSIN) || 
      (img -> code () == COU_EXPRPOW) || 
      (img -> code () == COU_EXPREXP) || 
      (img -> code () == COU_EXPRLOG) || 
      (img -> code () == COU_EXPRCOS)) {

    //    fprintf (stderr, " ==> "); w -> print (std::cerr); fprintf (stderr, "\n");

    expression *lbe, *ube;

    expression *indep = img -> Argument ();

    if (!indep) 
      indep = img -> ArgList () [0];//getFixVar ();

    int xi = indep -> Index ();

    //    fprintf (stderr, "looking into w_%d = f (x_%d)\n", 
    //       w -> Index (), xi);

    indep -> getBounds (lbe, ube);

    CouNumber lb = (*lbe) (),
      ub = (*ube) ();

    delete lbe;
    delete ube;

    if (xi >= 0) {

      CouNumber curx = cg -> Problem () -> X (xi);

#define N_STEPS 100

      // plot function

      if (first_draw) {

        first_draw = false;

        for (CouNumber x = lb; 
             x <= ub + COUENNE_EPS; 
             x += ((ub - lb) / N_STEPS)) {

          cg -> Problem () -> X () [xi] = x;
          
          CouNumber y = (*img) ();

          if (y > maxY) maxY = y;
          if (y < minY) minY = y;

          fprintf (stderr, "%.12e %.12e\n", x, y);
        }

        maxY += (maxY-minY) / 20;
        minY -= (maxY-minY) / 20;
      }
        
      lb -= (ub-lb) / 20;
      ub += (ub-lb) / 20;

      // plot lines defining constraint (only for cuts involving at
      // most two variables (that is, w is a unary function)
      for (int jj=j; jj < cs.sizeRowCuts (); jj++) {

        CouNumber lb0 = lb, 
          ub0 = ub;

        const double *el  = cs.rowCutPtr (jj) -> row (). getElements ();
        double  rhs = cs.rowCutPtr (jj) -> rhs ();

        if (fabs (el [1]) > COUENNE_EPS) {
          lb0 = CoinMax (lb, CoinMin ((rhs - el[0] * minY) / el [1], (rhs - el[0] * maxY) / el [1]));
          ub0 = CoinMin (ub, CoinMax ((rhs - el[0] * minY) / el [1], (rhs - el[0] * maxY) / el [1]));
        }

        fprintf (stderr, "#m=2,S=%d\n", (cs.rowCutPtr (jj) -> sense () == 'L') ? 10:11);

        fprintf (stderr, "%.12e %.12e\n", lb0, (rhs - el [1] * lb0) / el [0]);
        fprintf (stderr, "%.12e %.12e\n", ub0, (rhs - el [1] * ub0) / el [0]);
      }

      cg -> Problem () -> X () [xi] = curx;
      exit(0);
    }
  }
}