/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*    Copyright (C) 2002-2008 Konrad-Zuse-Zentrum                            */
/*                            fuer Informationstechnik Berlin                */
/*                                                                           */
/*  SCIP is distributed under the terms of the ZIB Academic License.         */
/*                                                                           */
/*  You should have received a copy of the ZIB Academic License.             */
/*  along with SCIP; see the file COPYING. If not email to scip@zib.de.      */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file   ConsNosubtour.cpp
 * @brief  C++ nosubtour constraint handler for the symmetric TSP
 * @author Timo Berthold
 * @author Kati Wolter
 */

/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/

#include <cassert>
#include <string>
#include <iostream>
#include "ConsNosubtour.h"
#include "objscip/objscip.h"

using namespace tsp;
using namespace scip;
using namespace std;

struct SCIP_ConsData
{
   GRAPH* graph;                          /**< the underlying graph */
};


/** checks given nosubtour constraints for feasibility of given solution or current solution */
static
void checkCons( 
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONS**        conss,              /**< array of constraints to process */
   int                nconss,             /**< number of constraints to process */
   SCIP_SOL*          sol,                /**< proposed solution, NULL for current LP or pseudo solution */
   SCIP_RESULT*       result              /**< pointer to store the result of the enforcing/checking call */
   )
{  
   assert(result != NULL);
   *result = SCIP_FEASIBLE;

   for( int i = 0; i < nconss; ++i )
   {
      SCIP_CONSDATA* consdata;
      GRAPH* graph;
      GRAPHNODE* node;
      GRAPHNODE* startnode;
      GRAPHEDGE* lastedge;
      GRAPHEDGE* edge;
      int tourlength;

      /* get underlying graph from constraint data structure */
      consdata = SCIPconsGetData(conss[i]);
      assert(consdata != NULL);
      graph = consdata->graph;
      assert(graph != NULL);   

      if(graph->nnodes <= 1)
         continue;

      startnode = &graph->nodes[0];
      tourlength = 0;
      lastedge = NULL;
      node = startnode;

      // follow the (sub?)tour until you come back to the startnode
      do
      {
         edge = node->first_edge;

         // look for an outgoing edge to proceed
         while( SCIPgetSolVal(scip, sol, edge->var) < 0.5 || edge->back == lastedge)
         {
            edge = edge->next;        
            assert(edge != NULL);
         }
         
         tourlength++;   
         node = edge->adjac;
         lastedge = edge;
      }
      while( node != startnode );

      assert(tourlength <= graph->nnodes);

      if( tourlength < graph->nnodes )
      {
         *result = SCIP_INFEASIBLE;
         return;
      }
   }
}

/** separates given nosubtour constraints: adds subtour elemination inequalities, 
 *  if violated by of given solution or current LP solution 
 */
static
SCIP_RETCODE sepaCons(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONS**        conss,              /**< array of constraints to process */
   int                nconss,             /**< number of constraints to process */
   SCIP_SOL*          sol,                /**< primal solution that should be separated */
   SCIP_RESULT*       result              /**< pointer to store the result of the separation call */
   )
{
   assert(result != NULL);

   *result = SCIP_DIDNOTFIND;

   for( int c = 0; c < nconss; ++c )
   {
      SCIP_CONSDATA* consdata;
      GRAPH* graph;
      SCIP_Bool** cuts;
      SCIP_Real cap;
      int ncuts;
      SCIP_Bool success;

      /* get underlying graph from constraint data structure */
      consdata = SCIPconsGetData(conss[c]);
      assert(consdata != NULL);
      graph = consdata->graph;
      assert(graph != NULL);
          
      // store the suggested, but infeasible solution into the capacity of the edges
      for( int i = 0; i < graph->nedges; i++)
      {
         cap = SCIPgetSolVal(scip, sol, graph->edges[i].var);
         graph->edges[i].rcap = cap;
         graph->edges[i].cap = cap;
         graph->edges[i].back->rcap = cap;
         graph->edges[i].back->cap = cap;   
      }
           
      SCIP_CALL( SCIPallocBufferArray(scip, &cuts, graph->nnodes) );
      for( int i = 0; i < graph->nnodes; i++)
      {
         SCIP_CALL( SCIPallocBufferArray(scip, &cuts[i], graph->nnodes) );
      }

      // try to find cuts
      success =  ghc_tree( graph, cuts, &ncuts, SCIPfeastol(scip) );
      if( success )
      { 
         int i = 0;

         // create a new cutting plane for every suitable arc (representing a cut with value < 2) of the Gomory Hu Tree
         while( i < ncuts )
         {
            SCIP_ROW* row; 
            SCIP_CALL( SCIPcreateEmptyRow(scip, &row, "subtour_elem", 2.0, SCIPinfinity(scip), FALSE, FALSE, TRUE) ); 

            SCIP_CALL( SCIPcacheRowExtensions(scip, row) );

            for( int j = 0; j < graph->nnodes; j++)
            { 
               // in gmincut the graph has been partitioned into two parts, represented by bools
               if( cuts[i][j] )
               {
                  GRAPHEDGE* edge = graph->nodes[j].first_edge;
                                        
                  // take every edge with nodes in different parts into account
                  while( edge != NULL )
                  {
                     if( !cuts[i][edge->adjac->id] )
                     {
                        SCIP_CALL( SCIPaddVarToRow(scip, row, edge->var, 1.0) ); 
                     }
                     edge = edge->next;
                  }
               }
            } 

            SCIP_CALL( SCIPflushRowExtensions(scip, row) );

            // add cut
            if( SCIPisCutEfficacious(scip, sol, row) )
            {
               SCIP_CALL( SCIPaddCut(scip, sol, row, FALSE) );
               *result = SCIP_SEPARATED;    
            }
            SCIP_CALL( SCIPreleaseRow(scip, &row) );
               
            i++;
         }            
      }

      for( int i = graph->nnodes - 1; i >= 0; i-- )
         SCIPfreeBufferArray( scip, &cuts[i] );
      SCIPfreeBufferArray( scip, &cuts );                
   }

   return SCIP_OKAY;
}


/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/


/** feasibility check method of constraint handler for primal solutions */
SCIP_RETCODE ConsNosubtour::scip_check(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS**        conss,              /**< array of constraints to process */
   int                nconss,             /**< number of constraints to process */
   SCIP_SOL*          sol,                /**< the solution to check feasibility for */
   SCIP_Bool          checkintegrality,   /**< has integrality to be checked? */
   SCIP_Bool          checklprows,        /**< have current LP rows to be checked? */
   SCIP_Bool          printreason,        /**< should the reason for the violation be printed? */
   SCIP_RESULT*       result              /**< pointer to store the result of the feasibility checking call */
   )
{
   *result = SCIP_INFEASIBLE;

   /* Nosubtour constraints are checked  after the node degree and integrality constraints.
    * If one of them was infeasible, we would not have been called
    */
   checkCons(scip, conss, nconss, sol, result);
         
   return SCIP_OKAY;
}

/** constraint enforcing method of constraint handler for LP solutions */
SCIP_RETCODE ConsNosubtour::scip_enfolp(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS**        conss,              /**< array of constraints to process */
   int                nconss,             /**< number of constraints to process */
   int                nusefulconss,       /**< number of useful (non-obsolete) constraints to process */
   SCIP_Bool          solinfeasible,      /**< was the solution already declared infeasible by a constraint handler? */
   SCIP_RESULT*       result              /**< pointer to store the result of the enforcing call */
   )
{
   *result = SCIP_INFEASIBLE;

   /* Nosubtour constraints are enforced after the node degree and integrality constraints.
    * We will only check feasibility, when the others are satisfied. 
    */
   if( !solinfeasible )
      checkCons(scip, conss, nconss, NULL, result);
         
   return SCIP_OKAY;
}

/** constraint enforcing method of constraint handler for pseudo solutions */
SCIP_RETCODE ConsNosubtour::scip_enfops(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS**        conss,              /**< array of constraints to process */
   int                nconss,             /**< number of constraints to process */
   int                nusefulconss,       /**< number of useful (non-obsolete) constraints to process */
   SCIP_Bool          solinfeasible,      /**< was the solution already declared infeasible by a constraint handler? */
   SCIP_Bool          objinfeasible,      /**< is the solution infeasible anyway due to violating lower objective bound? */
   SCIP_RESULT*       result              /**< pointer to store the result of the enforcing call */
   )
{
   *result = SCIP_INFEASIBLE;

   /* Nosubtour constraints are enforced after the node degree and integrality constraints.
    * We will only check feasibility, when the others are satisfied.
    */
   if( !solinfeasible )
      checkCons(scip, conss, nconss, NULL, result);
         
   return SCIP_OKAY;
}

/** separation method of constraint handler for LP solution */
SCIP_RETCODE ConsNosubtour::scip_sepalp(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS**        conss,              /**< array of constraints to process */
   int                nconss,             /**< number of constraints to process */
   int                nusefulconss,       /**< number of useful (non-obsolete) constraints to process */
   SCIP_RESULT*       result              /**< pointer to store the result of the separation call */
   )
{
   *result = SCIP_DIDNOTFIND;
   
   SCIP_CALL( sepaCons(scip, conss, nconss, NULL, result) );

   return SCIP_OKAY;
}

/** separation method of constraint handler for arbitrary primal solution */
SCIP_RETCODE ConsNosubtour::scip_sepasol(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS**        conss,              /**< array of constraints to process */
   int                nconss,             /**< number of constraints to process */
   int                nusefulconss,       /**< number of useful (non-obsolete) constraints to process */
   SCIP_SOL*          sol,                /**< primal solution that should be separated */
   SCIP_RESULT*       result              /**< pointer to store the result of the separation call */
   )
{
   *result = SCIP_DIDNOTFIND;

   SCIP_CALL( sepaCons(scip, conss, nconss, sol, result) );

   return SCIP_OKAY;
}


/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/


/** frees specific constraint data */
SCIP_RETCODE ConsNosubtour::scip_delete(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS*         cons,               /**< the constraint belonging to the constraint data */
   SCIP_CONSDATA**    consdata            /**< pointer to the constraint data to free */
   )
{
   assert(consdata != NULL);

   release_graph(&(*consdata)->graph);
   SCIPfreeMemory(scip, consdata);

   return SCIP_OKAY;
}


/** transforms constraint data into data belonging to the transformed problem */
SCIP_RETCODE ConsNosubtour::scip_trans(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS*         sourcecons,         /**< source constraint to transform */
   SCIP_CONS**        targetcons          /**< pointer to store created target constraint */
   )
{
   SCIP_CONSDATA* sourcedata;
   SCIP_CONSDATA* targetdata;

   sourcedata = SCIPconsGetData(sourcecons);

   SCIP_CALL( SCIPallocMemory(scip,&targetdata) );
   targetdata->graph = sourcedata->graph;
   capture_graph(targetdata->graph);

   /* create target constraint */
   SCIP_CALL( SCIPcreateCons(scip, targetcons, SCIPconsGetName(sourcecons), conshdlr, targetdata,
         SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons), SCIPconsIsEnforced(sourcecons),
         SCIPconsIsChecked(sourcecons), SCIPconsIsPropagated(sourcecons),  SCIPconsIsLocal(sourcecons),
         SCIPconsIsModifiable(sourcecons), SCIPconsIsDynamic(sourcecons), SCIPconsIsRemovable(sourcecons),
         SCIPconsIsStickingAtNode(sourcecons)) );

   return SCIP_OKAY;
}

/** variable rounding lock method of constraint handler */
SCIP_RETCODE ConsNosubtour::scip_lock(
   SCIP*              scip,               /**< SCIP data structure */
   SCIP_CONSHDLR*     conshdlr,           /**< the constraint handler itself */
   SCIP_CONS*         cons,               /**< the constraint that should lock rounding of its variables, or NULL if the
                                           *   constraint handler does not need constraints */
   int                nlockspos,          /**< no. of times, the roundings should be locked for the constraint */
   int                nlocksneg           /**< no. of times, the roundings should be locked for the constraint's negation */
   )
{
   SCIP_CONSDATA* consdata;
   GRAPH* g;
   
   consdata = SCIPconsGetData(cons);
   assert(consdata != NULL);
   
   g = consdata->graph;
   assert(g != NULL);

   for( int i = 0; i < g->nedges; ++i )
   {
      SCIP_CALL( SCIPaddVarLocks(scip, g->edges[i].var, nlocksneg, nlockspos) );
   }

   return SCIP_OKAY;
}


/** creates and captures a nosubtour constraint */
SCIP_RETCODE tsp::SCIPcreateConsNosubtour(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS**           cons,               /**< pointer to hold the created constraint */
   const char*           name,               /**< name of constraint */
   GRAPH*                graph,              /**< the underlying graph */
   SCIP_Bool             initial,            /**< should the LP relaxation of constraint be in the initial LP? */
   SCIP_Bool             separate,           /**< should the constraint be separated during LP processing? */
   SCIP_Bool             enforce,            /**< should the constraint be enforced during node processing? */
   SCIP_Bool             check,              /**< should the constraint be checked for feasibility? */
   SCIP_Bool             propagate,          /**< should the constraint be propagated during node processing? */
   SCIP_Bool             local,              /**< is constraint only valid locally? */
   SCIP_Bool             modifiable,         /**< is constraint modifiable (subject to column generation)? */
   SCIP_Bool             dynamic,            /**< is constraint dynamic? */
   SCIP_Bool             removable           /**< should the constraint be removed from the LP due to aging or cleanup? */
   )
{
   SCIP_CONSHDLR* conshdlr;
   SCIP_CONSDATA* consdata;

   /* find the subtour constraint handler */
   conshdlr = SCIPfindConshdlr(scip, "nosubtour");
   if( conshdlr == NULL )
   {
      SCIPerrorMessage("nosubtour constraint handler not found\n");
      return SCIP_PLUGINNOTFOUND;
   }

   /* create constraint data */
   SCIP_CALL( SCIPallocMemory( scip, &consdata) );
   consdata->graph = graph;
   capture_graph( consdata->graph );

   /* create constraint */
   SCIP_CALL( SCIPcreateCons(scip, cons, name, conshdlr, consdata, initial, separate, enforce, check, propagate,
         local, modifiable, dynamic, removable, FALSE) );

   return SCIP_OKAY;
}