CouenneFeasPump.cpp

Go to the documentation of this file.

/* $Id: CouenneFeasPump.cpp 1073 2014-04-11 01:38:17Z pbelotti $
 *
 * Name:    CouenneFeasPump.cpp
 * Authors: Pietro Belotti
 *          Timo Berthold, ZIB Berlin
 * Purpose: Implement the Feasibility Pump heuristic class
 * Created: August 5, 2009
 * 
 * This file is licensed under the Eclipse Public License (EPL)
 */

#include "CbcModel.hpp"
#include "CoinTime.hpp"
#include "CoinHelperFunctions.hpp"

#include "CouenneExprAux.hpp"
#include "CouenneFeasPump.hpp"
#include "CouenneProblem.hpp"
#include "CouenneProblemElem.hpp"
#include "CouenneCutGenerator.hpp"
#include "CouenneTNLP.hpp"
#include "CouenneFPpool.hpp"
#include "CouenneRecordBestSol.hpp"

#ifdef COIN_HAS_SCIP
/* general SCIP includes */
#include "scip/scip.h"
#include "scip/cons_linear.h"
#include "scip/scipdefplugins.h"
#endif

using namespace Ipopt;
using namespace Couenne;

void printDist   (CouenneProblem *p, const double *iSol, double *nSol);
void printCmpSol (CouenneProblem *p, const double *iSol, double *nSol, int direction);

// Solve
int CouenneFeasPump::solution (double &objVal, double *newSolution) {

  const int depth = (model_ && (model_ -> currentNode ())) ? model_ -> currentNode () -> depth () : 0;

  double time0 = CoinCpuTime();

  match_ = NULL;

  if ((nCalls_ == 0) ||                                   // check upper limit on number of runs
      //(problem_ -> nIntVars () <= 0) ||                 // feas pump on NLP? Why not?
      (CoinCpuTime () >= problem_ -> getMaxCpuTime ()) || // don't start if time is out
      ((numberSolvePerLevel_ >= 0) &&                     // stop FP after a certain level
       (CoinDrand48 () > 1. / CoinMax                     // decided randomly and inversely proportional
        (1., (double) ((CoinMax (0, depth - numberSolvePerLevel_ + 1)) *   // to BB tree depth
                                   (depth - numberSolvePerLevel_ + 1))))))
    return 0;

  if (nCalls_ > 0) // only decrease it if positive. If -1, no limit on # calls
    --nCalls_;

  problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "==================================================== FP: BEGIN\n");
  problem_ -> Jnlst () -> Printf (J_ERROR,   J_NLPHEURISTIC, "[FeasPump] Initializing\n");

  // Solve NLP once at the beginning ////////////////////////
  //
  // First, it possibly provides an NLP feasible point (the LP
  // provides a neither NLP nor integer point)
  //
  // Second, we need it to compile the Lagrangian Hessian used as a
  // distance

  if (!nlp_) // first call (in this run of FP). Create NLP
    nlp_ = new CouenneTNLP (problem_);

  //problem_ -> domain () -> push (*(problem_ -> domain () -> current ()));
  problem_ -> domain () -> push (model_ -> solver ());
  // NOTE: FP should only get bounds from the current BB node

  // Initial Bound Tightening: since Cbc prefers (reasonably) to call
  // a heuristic before any cut separator and because BT is instead
  // quite fast, let's see if we can reduce bounds to get a tighter
  // MILP relaxation.

  t_chg_bounds *chg_bds = new t_chg_bounds [problem_ -> nVars ()];

  for (int i=problem_ -> nVars (); i--;) {
    chg_bds [i]. setLower (t_chg_bounds::CHANGED);
    chg_bds [i]. setUpper (t_chg_bounds::CHANGED);
  }

  bool is_still_feas = problem_ -> btCore (chg_bds);

  delete [] chg_bds;

  // Now that bounds are possibly reduced, put current point within
  // bounds

  double
    *lb  = problem_ -> domain () -> lb (),
    *ub  = problem_ -> domain () -> ub (),
    *sol = problem_ -> domain () -> x  ();

  for (int i=problem_ -> nVars (); i--;)

    if      (sol [i] < lb [i]) sol [i] = lb [i];
    else if (sol [i] > ub [i]) sol [i] = ub [i];

  // fix integer coordinates of current (MINLP feasible!) solution
  // and set it as initial (obviously NLP feasible) solution

  nlp_ -> setInitSol (sol);


  if ((multHessNLP_  != 0.) || 
      (multHessMILP_ != 0.))
    nlp_ -> getSaveOptHessian () = true;

  if (problem_ -> Jnlst () -> ProduceOutput (J_ALL, J_NLPHEURISTIC)) {
    printf ("initial: solving NLP:---------\n");
    problem_ -> print ();
    printf ("---------------------\n");
  }

  problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: Initial NLP... "); fflush (stdout);

  // Solve with original objective function                           //       /|---------+
  ApplicationReturnStatus status = app_ -> OptimizeTNLP (nlp_);       //      < |   NLP   |
                                                                      //       \|---------+

  problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "done\n"); fflush (stdout);

  if (nlp_ -> getSolution () && (problem_ -> Jnlst () -> ProduceOutput (J_ALL, J_NLPHEURISTIC))) { // check if non-NULL
    printf ("######################## NLP solution (init):\n");
    for (int i=0; i< problem_ -> nVars ();) {
      printf ("%+e ", nlp_ -> getSolution () [i]);
      if (!(++i % 15)) printf ("\n");
    }
    if ((problem_ -> nVars () - 1) % 5) printf ("\n");
  }


  //problem_ -> domain () -> pop ();

  if ((status != Solve_Succeeded) && (status != Solved_To_Acceptable_Level))
    problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "Feasibility Pump: error in initial NLP problem\n");

  if ((multHessNLP_  != 0.) ||
      (multHessMILP_ != 0.))
    nlp_ -> getSaveOptHessian () = false;

  // This FP works as follows:
  //
  // obtain current NLP solution xN or, if none available, current LP solution
  //
  // repeat {
  //
  //   1) compute MILP solution(s) {xI} H_1-closest to xN
  //   2) insert them in pool
  //   3) consider the most promising one xI in the whole pool
  //
  //   if xI is MINLP-feasible (hack incumbent callback)
  //     update cutoff with min obj among MINLP-feasible
  //     run BT
  //   [else // a little out of the FP scheme
  //     apply linearization cuts to MILP]
  //
  //   select subset of variables (cont or integer), fix them
  //
  //   compute NLP solution xN that is H_2-closest to xI
  //
  //   if MINLP feasible
  //     update cutoff with min obj among MINLP-feasible
  //     run BT
  //
  // } until exit condition satisfied
  //
  // fix integer variables
  // resolve NLP

  CouNumber 
    *nSol = NULL, // solution of the nonlinear problem
    *iSol = NULL, // solution of the IP problem
    *best = NULL; // best solution found so far

  int
    niter    = 0,   // iteration counter
    nsuciter = 0,   // counter for consecutive successful applications of one MILP solving method
    retval   = 0,   // 1 if found a better solution
    nSep     = 0,   // If separation short circuit, max # of consecutive separations
    objInd   = problem_ -> Obj (0) -> Body () -> Index ();

  if (nlp_ -> getSolution ())
    nSol = CoinCopyOfArray (nlp_ -> getSolution (), problem_ -> nVars ());

  bool isNlpFeas = problem_ -> checkNLP0 (nSol, 
                                          nSol [objInd], 
                                          true, 
                                          false, // don't care about obj
                                          true,  // stop at first violation
                                          true); // checkAll

  //                      _                   _
  //                     (_)                 | |
  //  _ __ ___     __ _   _    _ __          | |    ___     ___    _ __ 
  // | '_ ` _ \   / _` | | |  | '_ \         | |   / _ \   / _ \  | '_ `.
  // | | | | | | | (_| | | |  | | | |        | |  | (_) | | (_) | | |_) |
  // |_| |_| |_|  \__,_| |_|  |_| |_|        |_|   \___/   \___/  | .__/
  //                                                              | |

  // copy bounding box and current solution to the problem (better
  // linearization cuts). Note that this push is pop()'d at the end
  // of the main routine

  // use new bounds found in BT earlier

  //problem_ -> domain () -> push (model_ -> solver ());

  expression *originalObjective = problem_ -> Obj (0) -> Body ();

  double
    save_mDN = multDistNLP_,
    save_mHN = multHessNLP_,
    save_mON = multObjFNLP_,
    save_mDM = multDistMILP_,
    save_mHM = multHessMILP_,
    save_mOM = multObjFMILP_;

  if (!isNlpFeas) do { // only do this if NLP solution was not also MINLP

    if (CoinCpuTime () > problem_ -> getMaxCpuTime ())
      break;

    if (niter) {
      multDistNLP_  *= fabs (save_mDN); // update within loop
      multHessNLP_  *= fabs (save_mHN);
      multObjFNLP_  *= fabs (save_mON);
      multDistMILP_ *= fabs (save_mDM);
      multHessMILP_ *= fabs (save_mHM);
      multObjFMILP_ *= fabs (save_mOM);
    }

    problem_ -> Jnlst () -> Printf (J_ERROR, J_NLPHEURISTIC, "[FeasPump] Iteration %d [%gs]\n", niter, CoinCpuTime() - time0);

    // INTEGER PART /////////////////////////////////////////////////////////

    // Solve IP using nSol as the initial point to minimize weighted
    // l-1 distance from. If nSol==NULL, the MILP is created using the
    // original milp's LP solution.

    bool bad_IP_sol = false;

    double z = solveMILP (nSol, iSol, niter, &nsuciter);

#if 0
    if (milp_ && milp_ -> getColLower ()) {
      printf ("fixed in milp_:    ");
      for (int i=0; i<problem_ -> nVars (); ++i)
        if (fabs (milp_ -> getColUpper () [i] - 
                  milp_ -> getColLower () [i]) < COUENNE_EPS)
          printf ("%d ", i);
    }

    problem_ -> domain () -> push (model_ -> solver ());
    printf ("\nfixed in model_:   ");
    for (int i=0; i<problem_ -> nVars (); ++i)
      if (fabs (problem_ -> Ub (i) - problem_ -> Lb (i)) < COUENNE_EPS)
        printf ("%d ", i);
    problem_ -> domain () -> pop ();

    problem_ -> domain () -> push (*(problem_ -> domain () -> current ()));
    printf ("\nfixed in problem_: ");
    for (int i=0; i<problem_ -> nVars (); ++i)
      if (fabs (problem_ -> Ub (i) - problem_ -> Lb (i)) < COUENNE_EPS)
        printf ("%d ", i);
    problem_ -> domain () -> pop ();
    printf ("\n");
#endif

    // if no MILP solution was found, bail out

    bool try_again = false;

    if (!iSol || z >= COIN_DBL_MAX/2) {

      problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: could not find IP solution\n");

      // find non-tabu solution in the solution pool
      while (!pool_ -> Set (). empty ()) {

        // EXTRACT the closest (to nSol) IP solution from the pool
        pool_ -> findClosestAndReplace (iSol, iSol, problem_ -> nVars ());

        CouenneFPsolution newSol (problem_, iSol);

        // we found a solution that is not in the tabu list
        if (tabuPool_ . find (newSol) == tabuPool_ . end ()) {

          //tabuPool_ . insert (newSol);

          try_again = true;
          break;
        }
      }

      if (!try_again) { // try moving around current solution 

        bad_IP_sol = true;

        // SCIP could not find a MILP solution and we're somewhat
        // locked. Round current solution and feed it to the NLP. This
        // is better than just bailing out.

        int n = problem_ -> nVars ();

        if (!iSol)
          iSol = new double [n];

        for (int i=0; i<n; i++)
          iSol [i] = (problem_ -> Var (i) -> isInteger ()) ? 
            COUENNE_round (nSol [i]) : 
            nSol [i];
      }
 
      // if (!try_again) { // nothing to do, bail out   
      //        problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: could not find from pool either, bailing out\n");
      //        break;
      // }
    }

    bool isChecked = false;

    // If a solution was found, but is in the tabu list, two choices:
    //
    // 1) the pool is empty: do a round of cuts and repeat;
    //
    // 2) the pool is nonempty: extract the best solution from the
    //    pool and use it instead

    problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: %d solutions in pool, %d in tabu list\n", pool_ -> Set (). size (), tabuPool_ . size ());

    CouenneFPsolution checkedSol (problem_, iSol, false); // false is for not allocating space for this

    if (tabuPool_. find (checkedSol) == tabuPool_ . end ())

      tabuPool_. insert (CouenneFPsolution (problem_, iSol)); // only insertion to tabu pool: we check its feasibility now

    else {

      problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: found solution is tabu\n");

      // Current solution was visited before. Replace it with another
      // MILP solution from the pool, if any.

      if         (tabuMgt_ == FP_TABU_NONE) break;

      else   if ((tabuMgt_ == FP_TABU_POOL) && !(pool_ -> Set (). empty ())) {

        // try to find non-tabu solution in the solution pool

        do {

          // retrieve the top solution from the pool
          pool_ -> findClosestAndReplace (iSol, nSol, problem_ -> nVars ());

          CouenneFPsolution newSol (problem_, iSol);

          // we found a solution that is not in the tabu list
          if (tabuPool_. find (newSol) == tabuPool_ . end ()) {

            tabuPool_ . insert (newSol);
            break;
          }

          // the pool is empty -> just round
          if (pool_ -> Set ().empty ())
            {
              int n = problem_ -> nVars ();

              if (!iSol)
                iSol = new double [n];

              for (int i=0; i<n; i++)
                iSol [i] = (problem_ -> Var (i) -> isInteger ()) ? 
                  COUENNE_round (nSol [i]) : 
                  nSol [i];

              // delete[] iSol;
              // iSol = NULL;
            }

        } while( !pool_ -> Set ().empty() );

      } else if (((tabuMgt_ == FP_TABU_CUT)   ||  
                  ((pool_ -> Set (). empty ()) && iSol))) {

        OsiCuts cs;

        problem_   -> domain () -> push (problem_ -> nVars (), iSol, milp_ -> getColLower (), milp_ -> getColUpper (), false); // don't copy vectors
        couenneCG_ -> genRowCuts (*milp_, cs, 0, NULL); // remaining three arguments NULL by default
        problem_   -> domain () -> pop ();

        if (cs.sizeRowCuts ()) {

          milp_ -> applyCuts (cs);

          if (nSep++ < nSepRounds_)
            continue;

        } else break; // nothing left to do, just bail out

      } else if ((tabuMgt_ == FP_TABU_PERTURB) && iSol) {

        // perturb solution     

        const CouNumber 
          *lb = milp_ -> getColLower (),
          *ub = milp_ -> getColUpper ();

        double 
          downMoves = 0.,
          upMoves   = 0.;

        int startIndex = (int) floor (CoinDrand48 () * problem_ -> nOrigVars ());

        for (int ii = problem_ -> nOrigVars (); ii--; lb++, ub++) {

          if (problem_ -> Var (ii) -> Multiplicity () <= 0)
            continue;

          // make perturbation start from pseudo-random points

          int i = (startIndex + ii) % problem_ -> nOrigVars ();

          if (problem_ -> Var (i) -> isInteger ()) {

            double
              rnd  = CoinDrand48 (), 
              down = 0.,
              up   = 1.;

            // if there is room on the left (resp. right) of the
            // solution, consider moving down (resp. up). Make moves
            // less likely as we don't want to perturb too many
            // variables

#define RND_DECR_EXPONENT .5

            if (iSol [i] >= lb [i] + 1.) down =      1. / pow (1. + (downMoves += 1.), RND_DECR_EXPONENT);
            if (iSol [i] <= ub [i] - 1.) up   = 1. - 1. / pow (1. + (upMoves   += 1.), RND_DECR_EXPONENT);

            if      (rnd < down) iSol [i] -= 1.;
            else if (rnd > up)   iSol [i] += 1.;
          }
        }
      }
    }

    if (!iSol)
      break;

    problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: checking IP solution for feasibility\n");

    isChecked = problem_ -> checkNLP0 (iSol, z, true, 
                                       false, // don't care about obj
                                       true,  // stop at first violation
                                       true); // checkAll

    // Possible improvement: if IP-infeasible or if z does not improve
    // the best solution found so far, then try the following:
    //
    // 1) Fix integer variables to IP solution's corresponding components
    // 2) Solve restriction with original obj
    // 3) While not found (MI)NLP solution or
    //          solution MINLP infeasible or
    //          z not better
    // 4)   Get solution x* from pool
    // 5)   Solve with original objective
    // 6) Done
    // 7) If found solution, set it, otherwise keep previously found (IP-feasible) one

    if ((!isChecked ||                     // not MINLP feasible
         z > problem_ -> getCutOff ())) {  // not improving

      problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: infeasible/non-improving (feas==%d, z=%g, cutoff=%g), looping on pool\n", isChecked, z, problem_->getCutOff ());

      bool try_again;

      do {

        try_again = false;

        if (CoinCpuTime () > problem_ -> getMaxCpuTime ())
          break;

        // Check if fixing integers to those in iSol yields an
        // infeasible problem. If so, don't optimize
        if (fixIntVariables (iSol)) {

          nlp_ -> setInitSol (iSol);

          if (problem_ -> Jnlst () -> ProduceOutput (J_ALL, J_NLPHEURISTIC)) {
            printf ("----------------------- Solving NLP:\n");
            problem_ -> print ();
            printf ("-----------------------\n");
          }

          status = app_ -> OptimizeTNLP (nlp_);

          if (nlp_ -> getSolution ()) { // check if non-NULL
            if  (nSol)  CoinCopyN       (nlp_ -> getSolution (), problem_ -> nVars (), nSol);
            else nSol = CoinCopyOfArray (nlp_ -> getSolution (), problem_ -> nVars ());
          }

          if (nlp_ -> getSolution () && (problem_ -> Jnlst () -> ProduceOutput (J_ALL, J_NLPHEURISTIC))) { // check if non-NULL
            printf ("######################## NLP solution (loop through pool):\n");
              for (int i=0; i< problem_ -> nVars ();) {
                printf ("%+e ", nlp_ -> getSolution () [i]);
                if (!(++i % 15)) printf ("\n");
              }
          }

          z = nlp_ -> getSolValue ();

          if (z < problem_ -> getCutOff ()) // don't waste time if not improving

            isChecked = problem_ -> checkNLP0 (nSol, z, true, 
                                               false, // don't care about obj
                                               true,  // stop at first violation
                                               true); // checkAll

          if (isChecked && (z < problem_ -> getCutOff ())) 
            break;
        } 

        // ONLY do the following if we didn't get a good IP solution
        // earlier (i.e. if we had to round)

        // find non-tabu solution in the solution pool
        if (bad_IP_sol)
          while (!pool_ -> Set (). empty ()) {

            // EXTRACT the closest (to nSol) IP solution from the pool
            pool_ -> findClosestAndReplace (iSol, nSol, problem_ -> nVars ());

            CouenneFPsolution newSol (problem_, iSol);

            // we found a solution that is not in the tabu list
            if (tabuPool_ . find (newSol) == tabuPool_ . end ()) {

              tabuPool_ . insert (newSol);
              try_again = true;
              break;
            }
          } 

      } while (try_again);
    }

    // Whatever the previous block yielded, we are now going to check
    // if we have a new solution, and if so we save it.

    if (isChecked) {

      problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: IP solution is MINLP feasible\n");

      // solution is MINLP feasible! Save it.

      retval = 1;
      objVal = z;

      // Found a MINLP-feasible solution, but to keep diversity do NOT
      // use the best available. Just use this.
      //
      // #ifdef FM_CHECKNLP2
      // #  ifdef FM_TRACE_OPTSOL
      //       problem_->getRecordBestSol()->update();
      //       best = problem_->getRecordBestSol()->getSol();
      //       objVal = problem_->getRecordBestSol()->getVal();
      // #  else /* not FM_TRACE_OPTSOL */
      //       best = problem_->getRecordBestSol()->getModSol(problem_->nVars());
      //       objVal = z;
      // #  endif /* not FM_TRACE_OPTSOL */
      // #else /* not FM_CHECKNLP2 */
      // #  ifdef FM_TRACE_OPTSOL
      //       problem_->getRecordBestSol()->update(iSol, problem_->nVars(),
      //                                           z, problem_->getFeasTol());
      //       best = problem_->getRecordBestSol()->getSol();
      //       objVal = problem_->getRecordBestSol()->getVal();
      // #  else /* not FM_TRACE_OPTSOL */

      best = iSol;

      //       objVal = z;
      // #  ifdef FM_TRACE_OPTSOL
      //       problem_->getRecordBestSol()->update();
      //       best = problem_->getRecordBestSol()->getSol();
      //       objVal = problem_->getRecordBestSol()->getVal();
      // #  endif /* not FM_TRACE_OPTSOL */
      // #endif /* not FM_CHECKNLP2 */

      if (z < problem_ -> getCutOff ()) {

        problem_ -> setCutOff (objVal);

        t_chg_bounds *chg_bds = NULL;

        if (objInd >= 0) {
          chg_bds = new t_chg_bounds [problem_ -> nVars ()];
          chg_bds [objInd].setUpper (t_chg_bounds::CHANGED); 
        }

        // If, with new cutoff, bound tightening clears the whole
        // feasible set, stop
        bool is_still_feas = problem_ -> btCore (chg_bds);

        if (chg_bds)
          delete [] chg_bds;

        // don't tighten MILP if BT says it's infeasible

        if (!is_still_feas)
          break;

        // Update lb/ub on milp and nlp here
        const CouNumber
          *plb = problem_ -> Lb (),
          *pub = problem_ -> Ub (),
          *mlb = milp_    -> getColLower (),
          *mub = milp_    -> getColUpper ();

        for (int i=problem_ -> nVars (), j=0; i--; ++j, ++plb, ++pub) {

          bool neglect = problem_ -> Var (j) -> Multiplicity () <= 0;

          // neglect'ed variables don't appear in the problem, so MILP can fix them to zero
          if (*plb > *mlb++) milp_ -> setColLower (j, neglect ? 0. : *plb);
          if (*pub < *mub++) milp_ -> setColUpper (j, neglect ? 0. : *pub);
        }
      }

      break;

    } else {

      problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: IP solution NOT MINLP feasible\n");

      if (milpCuttingPlane_ == FP_CUT_EXTERNAL || 
          milpCuttingPlane_ == FP_CUT_POST) {

        // Solution is IP- but not MINLP feasible: it might get cut by
        // linearization cuts. If so, add a round of cuts and repeat.

        OsiCuts cs;

        problem_   -> domain () -> push (milp_);
        couenneCG_ -> genRowCuts (*milp_, cs, 0, NULL); // remaining three arguments NULL by default
        problem_   -> domain () -> pop ();

        if (cs.sizeRowCuts ()) { 

          // the (integer, NLP infeasible) solution could be separated

          milp_ -> applyCuts (cs);

          // found linearization cut, now re-solve MILP (not quite a FP)
          if (milpCuttingPlane_ == FP_CUT_EXTERNAL && 
              nSep++ < nSepRounds_)
            continue;
        }
      }
    }

    //
    // reset number of separation during non-separating iteration
    //

    nSep = 0;

    // NONLINEAR PART ///////////////////////////////////////////////////////

    // fix some variables and solve the NLP to find a NLP (possibly
    // non-MIP) feasible solution

    z = solveNLP (iSol, nSol); 

    if ((nSol && iSol) &&
        (problem_ -> Jnlst () -> ProduceOutput (J_WARNING, J_NLPHEURISTIC))) {

      double dist = 0.;
      int nNonint = 0;

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

        exprVar *e = problem_ -> Var (i);

        if (e -> Multiplicity () <= 0)
          continue;

        if (e -> isInteger () &&
            (fabs (iSol [i] - ceil (iSol [i] - .5)) > COUENNE_EPS))
          ++nNonint;

        dist += 
          (iSol [i] - nSol [i]) * 
          (iSol [i] - nSol [i]);
      }

      printf ("FP: after NLP, distance %g, %d nonintegers\n", sqrt (dist), nNonint);
    }

    if (problem_ -> Jnlst () -> ProduceOutput (J_WARNING, J_NLPHEURISTIC)) {

      printDist   (problem_, iSol, nSol);
      printCmpSol (problem_, iSol, nSol, 0);
    }

    if (z > COIN_DBL_MAX/2) // something went wrong in the NLP, bail out
      break;

    // check if newly found NLP solution is also integer

    isChecked = false;

    if (nSol) {

      isChecked = problem_->checkNLP0 (nSol, z, true,
                                       false, // do not care about obj 
                                       true, // stopAtFirstViol
                                       true); // checkALL
    }

    if (nSol && isChecked) {

      problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: NLP solution is MINLP feasible\n");

      retval = 1;
      objVal = z;

#ifdef FM_TRACE_OPTSOL // - if ----------------------------
#ifdef FM_CHECKNLP2
      problem_->getRecordBestSol()->update();
#else
      problem_->getRecordBestSol()->update(nSol, problem_->nVars(), z, problem_->getFeasTol());
#endif
      best = problem_->getRecordBestSol()->getSol();
      objVal = problem_->getRecordBestSol()->getVal();
#else                  // - else --------------------------
#ifdef FM_CHECKNLP2
      best = problem_->getRecordBestSol()->getModSol(problem_ -> nVars ());
#else
      best   = nSol;
#endif
      objVal = z;
#endif                 // - endif -------------------------

      if (z < problem_ -> getCutOff ()) {

        problem_ -> setCutOff (objVal);

        t_chg_bounds *chg_bds = NULL;

        if (objInd >= 0) {
        
          chg_bds = new t_chg_bounds [problem_ -> nVars ()];
          chg_bds [objInd].setUpper (t_chg_bounds::CHANGED); 
        }

        // if bound tightening clears the whole feasible set, stop
        bool is_still_feas = problem_ -> btCore (chg_bds);

        if (chg_bds) 
          delete [] chg_bds;

        if (!is_still_feas)
          break;
      }
    }

    // update weights for hessian, distance, and objective, for both
    // MILP and NLP parts

    //multHessNLP_ *= save_multHessNLP_; // exponential reduction

  } while ((niter++ < maxIter_) && 
           (retval == 0));

  problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: out of the FP loop (feas=%d\n)", retval);

  // OUT OF THE LOOP ////////////////////////////////////////////////////////

  // If MINLP solution found,
  //
  // 1) restore original objective 
  // 2) fix integer variables
  // 3) resolve NLP

  if (nlp_)
    nlp_ -> setObjective (originalObjective);

  if (retval > 0) {

    problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: fix integer variables and rerun NLP\n");

    if (!nlp_) // first call (in this run of FP). Create NLP
      nlp_ = new CouenneTNLP (problem_);

    problem_ -> domain () -> push (*(problem_ -> domain () -> current ()));

    // fix integer coordinates of current (MINLP feasible!) solution
    // and, if feasible, set it as initial (obviously NLP feasible)
    // solution

    if (fixIntVariables (best)) {

      nlp_ -> setInitSol (best);


      //app_ -> Options () -> SetStringValue ("fixed_variable_treatment", "make_parameter");

      // Solve with original objective function
      status = app_ -> OptimizeTNLP (nlp_);


      if ((status != Solve_Succeeded) && 
          (status != Solved_To_Acceptable_Level))
 
        problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, 
                                        "Feasibility Pump: error in final NLP problem (due to fixing integer variables?)\n");

      if (nlp_ -> getSolution () && (problem_ -> Jnlst () -> ProduceOutput (J_ALL, J_NLPHEURISTIC))) { // check if non-NULL
        printf ("######################## NLP solution (post):\n");
          for (int i=0; i< problem_ -> nVars ();) {
            printf ("%+e ", nlp_ -> getSolution () [i]);
            if (!(++i % 15)) printf ("\n");
          }
      }

      // if found a solution with the last NLP, check & save it

      double z = nlp_ -> getSolValue ();

      // check if newly found NLP solution is also integer (unlikely...)
      bool isChecked = false;

      if (nSol) {

        problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: found nlp solution, check it\n");

        isChecked = problem_ -> checkNLP0 (nSol, z, true,
                                           false,
                                           true,
                                           true);
      }

      if (nSol &&
          isChecked &&
          (z < problem_ -> getCutOff ())) {

        problem_ -> Jnlst () -> Printf (J_WARNING, J_NLPHEURISTIC, "FP: feasible solution is improving\n");

#ifdef FM_TRACE_OPTSOL

#ifdef FM_CHECKNLP2
        problem_->getRecordBestSol()->update();
#else
        problem_->getRecordBestSol()->update(nSol, problem_->nVars(), z, problem_->getFeasTol());
#endif
        best = problem_->getRecordBestSol()->getSol();
        objVal = problem_->getRecordBestSol()->getVal();
#else

#ifdef FM_CHECKNLP2
        best = problem_->getRecordBestSol()->getModSol(problem_ -> nVars ());
#else
        best   = nSol;
#endif
        objVal = z;
#endif

        problem_ -> setCutOff (objVal);
      }
    }

    problem_ -> domain () -> pop ();

    if (problem_ -> Jnlst () -> ProduceOutput (J_WARNING, J_NLPHEURISTIC)) {
      printf ("FP: returning MINLP feasible solution:\n");
      printDist (problem_, best, nSol);
    }

    CoinCopyN (best, problem_ -> nVars (), newSolution);
  }

  // post-call fading update of coefficients

  multDistNLP_  = save_mDN * fadeMult_;
  multHessNLP_  = save_mHN * fadeMult_;
  multObjFNLP_  = save_mON * fadeMult_;
  multDistMILP_ = save_mDM * fadeMult_;
  multHessMILP_ = save_mHM * fadeMult_;
  multObjFMILP_ = save_mOM * fadeMult_;

  if (iSol) delete [] iSol;
  if (nSol) delete [] nSol;

  // release bounding box
  problem_ -> domain () -> pop ();

  // deleted at every call from Cbc, since it changes not only in
  // terms of variable bounds but also in terms of linearization cuts
  // added

  delete milp_;
  delete postlp_;
  milp_ = postlp_ = NULL;

  if (match_) {
    delete [] match_;
    match_ = NULL;
  }

  problem_ -> Jnlst () -> Printf 
    (J_WARNING, J_NLPHEURISTIC, "FP: done ===================\n");

  if (!retval)
    problem_ -> Jnlst () -> Printf 
      (J_ERROR, J_NLPHEURISTIC, "FP: No solution found\n");

  if (retval && (-2 == nCalls_)) {
    problem_ -> bonBase () -> options () -> SetNumericValue ("time_limit", 0.001, "couenne.");
    problem_ -> bonBase () -> setDoubleParameter (Bonmin::BabSetupBase::MaxTime, 0.001);
  }

  return retval;
}

// gather data on single solution for later printout
void compDistSingle (CouenneProblem *p,
                     int n,
                     const double *v,
                     double &norm,
                     int &nInfI,
                     int &nInfN,
                     double &infI,
                     double &infN) {

  p -> domain () -> push (n, v, NULL, NULL); // don't care about bounds as not used

  norm = infI = infN = 0.;
  nInfI = nInfN = 0;

  while (n--) {
    norm += (*v * *v);
    ++v;
  }
  v -= n;

  norm = sqrt (norm);

  for (std::vector <exprVar *>::iterator i = p -> Variables (). begin (); 
       i != p -> Variables (). end ();
       ++i) {

    if ((*i) -> Multiplicity () <= 0)
      continue;

    CouNumber 
      vval = (**i) ();

    if ((*i) -> isInteger ()) {

      double inf = CoinMax (vval - floor (vval + COUENNE_EPS),
                            ceil (vval - COUENNE_EPS) - vval);

      if (inf > COUENNE_EPS) {
        ++nInfI;
        infI += inf;
      }
    }

    if ((*i) -> Type () == AUX) {

      double 
        diff = 0.,
        fval = (*((*i) -> Image ())) ();

      if      ((*i) -> sign () != expression::AUX_GEQ) diff = CoinMax (diff, vval - fval);
      else if ((*i) -> sign () != expression::AUX_LEQ) diff = CoinMax (diff, fval - vval);

      if (diff > COUENNE_EPS) {
        ++nInfN;
        infN += diff;
      }
    }
  }

  p -> domain () -> pop ();
}


// print solutions and distances
void printDist (CouenneProblem *p, const double *iSol, double *nSol) {

  int nInfII = -1, nInfNI = -1, nInfIN = -1, nInfNN = -1;

  double 
    dist  = -1., 
    normI = -1., normN = -1., 
    infII = -1., infNI = -1.,
    infIN = -1., infNN = -1.;

  if (iSol) compDistSingle (p, p -> nVars (), iSol, normI, nInfII, nInfNI, infII, infNI);
  if (nSol) compDistSingle (p, p -> nVars (), nSol, normN, nInfIN, nInfNN, infIN, infNN);
  
  if (iSol && nSol) {

    dist = 0.;

    for (int i = p -> nVars (); i--;) {

      if (p -> Var (i) -> Multiplicity () <= 0)
        continue;

      dist += 
        (iSol [i] - nSol [i]) * 
        (iSol [i] - nSol [i]);
    }

    dist = sqrt (dist);
  }

  printf ("FP: ");
  printf ("IP norm i:%e n:%e dist %e #inf i:%4d n:%4d max inf i:%e n:%e ", normI, normN, dist, nInfII, nInfNI, infII, infNI);
  printf ("NL #inf i:%4d n:%4d max inf i:%e n:%e\n", nInfIN, nInfNN, infIN, infNN);
}


#define WRAP 3

void printCmpSol (CouenneProblem *p, const double *iSol, double *nSol, int direction) {

  int n = p -> nVars ();

  printf ("i:%p n:%p\nFP # ", 
          (void *) iSol, (void *) nSol);

  double 
    distance = 0.,
    diff;

  char c =
    direction < 0 ? '<' : 
    direction > 0 ? '>' : '-';

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

    if (p -> Var (i) -> Multiplicity () <= 0)
      continue;

    if (i && !(i % WRAP))
      printf ("\nFP # ");

    double 
      iS = iSol ? iSol [i] : 12345., // flag values, to be seen if something is wrong
      nS = nSol ? nSol [i] : 54321.;

    printf ("[%4d %+e -%c- %+e (%e)] ", 
            i, iS, c, nS,
            (iSol && nSol) ? fabs (iS - nS) : 0.);

    if (iSol && nSol) {
      diff = iS - nS;
      distance += (diff*diff);
    }
  }

  if (iSol && nSol) {

    distance = sqrt (distance);
    printf ("\n### distance: %e\n", distance);
  }
}