/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program and library             */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/*    Copyright (C) 2002-2022 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 visit scipopt.org.         */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/* uncomment for debug output: */
/* #define SCIP_DEBUG */

/**@file   cons_lop.c
 * @brief  constraint handler for linear ordering constraints
 * @author Marc Pfetsch
 *
 * We handle the following system of linear constraints:
 * - \f$ x_{ij} + x_{ji} = 1 \f$ for \f$i < j\f$                               (symmetry equations - added initially)
 * - \f$ x_{ij} + x_{jk} + x_{ki} \leq 2 \f$ for \f$i < j, i < k, j \neq k\f$  (triangle inequalities - separated)
 */

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

#include "cons_lop.h"

#include <assert.h>
#include <string.h>


/* constraint handler properties */
#define CONSHDLR_NAME          "lop"
#define CONSHDLR_DESC          "linear ordering constraint handler"
#define CONSHDLR_SEPAPRIORITY       100 /**< priority of the constraint handler for separation */
#define CONSHDLR_ENFOPRIORITY      -100 /**< priority of the constraint handler for constraint enforcing */
#define CONSHDLR_CHECKPRIORITY     -100 /**< priority of the constraint handler for checking feasibility */
#define CONSHDLR_SEPAFREQ             1 /**< frequency for separating cuts; zero means to separate only in the root node */
#define CONSHDLR_PROPFREQ             1 /**< frequency for propagating domains; zero means only preprocessing propagation */
#define CONSHDLR_EAGERFREQ          100 /**< frequency for using all instead of only the useful constraints in separation,
                                         *   propagation and enforcement, -1 for no eager evaluations, 0 for first only */
#define CONSHDLR_DELAYSEPA        FALSE /**< should separation method be delayed, if other separators found cuts? */
#define CONSHDLR_DELAYPROP        FALSE /**< should propagation method be delayed, if other propagators found reductions? */
#define CONSHDLR_NEEDSCONS         TRUE /**< should the constraint handler be skipped, if no constraints are available? */

#define CONSHDLR_PROP_TIMING       SCIP_PROPTIMING_BEFORELP


/** constraint data for linear ordering constraints */
struct SCIP_ConsData
{
   int                   n;                  /**< number of elements */
   SCIP_VAR***           vars;               /**< variables */
};


/** separate symmetry equations and triangle inequalities */
static
SCIP_RETCODE LOPseparate(
   SCIP*                 scip,               /**< SCIP pointer */
   SCIP_CONSHDLR*        conshdlr,           /**< constraint handler */
   int                   n,                  /**< number of elements */
   SCIP_VAR***           vars,               /**< n x n matrix of variables */
   SCIP_SOL*             sol,                /**< solution to be separated */
   int*                  nGen,               /**< output: pointer to store number of added rows */
   SCIP_Bool*            cutoff              /**< output: pointer to store whether we detected a cutoff */
   )
{
   char s[SCIP_MAXSTRLEN];
   int i;
   int j;
   int k;

   assert( scip != NULL );
   assert( vars != NULL );
   assert( nGen != NULL );
   assert( cutoff != NULL );

   *cutoff = FALSE;
   for (i = 0; i < n && ! (*cutoff); ++i)
   {
      for (j = i+1; j < n && ! (*cutoff); ++j)
      {
	 SCIP_Real valIJ;

	 valIJ = SCIPgetSolVal(scip, sol, vars[i][j]);

	 /* if symmetry equations are violated - should not be the case, if they are added in the beginning */
	 if ( ! SCIPisFeasEQ(scip, valIJ + SCIPgetSolVal(scip, sol, vars[j][i]), 1.0) )
	 {
	    SCIP_ROW *row;

	    (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "sym#%d#%d", i, j);

	    SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, s, 1.0, 1.0, FALSE, FALSE, TRUE) );
	    SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
	    SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i][j], 1.0) );
	    SCIP_CALL( SCIPaddVarToRow(scip, row, vars[j][i], 1.0) );
	    SCIP_CALL( SCIPflushRowExtensions(scip, row) );
#ifdef SCIP_DEBUG
	    SCIPdebug( SCIPprintRow(scip, row, NULL) );
#endif
	    SCIP_CALL( SCIPaddRow(scip, row, FALSE, cutoff) );
	    SCIP_CALL( SCIPreleaseRow(scip, &row));
	    ++(*nGen);

            if ( *cutoff )
               break;
	 }

	 /* check triangle inequalities */
	 for (k = i+1; k < n; ++k)
	 {
	    SCIP_Real sum;

	    if ( k == j )
	       continue;

	    sum = valIJ + SCIPgetSolVal(scip, sol, vars[j][k]) + SCIPgetSolVal(scip, sol, vars[k][i]);

	    /* if sum - 2.0 > 0, i.e., the cut is violated */
	    if ( SCIPisEfficacious(scip, sum - 2.0) )
	    {
	       SCIP_ROW *row;

	       (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "triangle#%d#%d#%d", i, j, k);

	       SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, s, -SCIPinfinity(scip), 2.0, FALSE, FALSE, TRUE) );
	       SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
	       SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i][j], 1.0) );
	       SCIP_CALL( SCIPaddVarToRow(scip, row, vars[j][k], 1.0) );
	       SCIP_CALL( SCIPaddVarToRow(scip, row, vars[k][i], 1.0) );
	       SCIP_CALL( SCIPflushRowExtensions(scip, row) );
#ifdef SCIP_DEBUG
	       SCIPdebug( SCIPprintRow(scip, row, NULL) );
#endif
	       SCIP_CALL( SCIPaddRow(scip, row, FALSE, cutoff) );
	       SCIP_CALL( SCIPreleaseRow(scip, &row));
	       ++(*nGen);

               if ( *cutoff )
                  break;
	    }
	 }
      }
   }

   return SCIP_OKAY;
}


/** copy method for constraint handler plugins (called when SCIP copies plugins) */
static
SCIP_DECL_CONSHDLRCOPY(conshdlrCopyLOP)
{  /*lint --e{715}*/
   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( valid != NULL );

   /* call inclusion method of constraint handler */
   SCIP_CALL( SCIPincludeConshdlrLOP(scip) );

   *valid = TRUE;

   return SCIP_OKAY;
}

/** frees specific constraint data */
static
SCIP_DECL_CONSDELETE(consDeleteLOP)
{  /*lint --e{715}*/
   int i;
   int n;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( cons != NULL );
   assert( consdata != NULL);
   assert( *consdata != NULL);
   assert( (*consdata)->vars != NULL );

   SCIPdebugMsg(scip, "deleting linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

   n = (*consdata)->n;
   for (i = 0; i < n; ++i)
      SCIPfreeBlockMemoryArray(scip, &((*consdata)->vars[i]), n); /*lint !e866*/
   SCIPfreeBlockMemoryArray(scip, &((*consdata)->vars), n);
   SCIPfreeBlockMemory(scip, consdata);

   return SCIP_OKAY;
}

/** deinitialization method of constraint handler (called before transformed problem is freed)
 *
 *  We output the final linear ordering.
 */
static
SCIP_DECL_CONSEXIT(consExitLOP)
{  /*lint --e{715}*/
   SCIP_SOL* sol;
   int c;
   int i;
   int j;
   int n;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );

   SCIPdebugMsg(scip, "exiting linear ordering constraint handler <%s>.\n", SCIPconshdlrGetName(conshdlr));

   /* avoid output for subscips */
   if ( SCIPgetSubscipDepth(scip) > 0 )
      return SCIP_OKAY;

   /* get best solution */
   sol = SCIPgetBestSol(scip);
   if ( sol == NULL )
      return SCIP_OKAY;

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_VAR*** vars;
      int* outdeg;
      int* indices;

      assert( conss != NULL );
      assert( conss[c] != NULL );
      SCIPdebugMsg(scip, "solution for for linear ordering constraint <%s>.\n", SCIPconsGetName(conss[c]));

      consdata = SCIPconsGetData(conss[c]);
      assert( consdata != NULL );
      assert( consdata->vars != NULL );
      n = consdata->n;
      vars = consdata->vars;

      SCIP_CALL( SCIPallocBufferArray(scip, &outdeg, n) );
      SCIP_CALL( SCIPallocBufferArray(scip, &indices, n) );

      /* compute out-degree */
      for (i = 0; i < n; ++i)
      {
	 int deg = 0;
	 for (j = 0; j < n; ++j)
	 {
	    SCIP_Real val;

	    if (j == i)
	       continue;

	    val = SCIPgetSolVal(scip, sol, vars[i][j]);
	    assert( SCIPisFeasIntegral(scip, val) );
	    if ( val < 0.5 )
	       ++deg;
	 }
	 outdeg[i] = deg;
	 indices[i] = i;
      }

      /* sort such that degrees are non-decreasing */
      SCIPsortIntInt(outdeg, indices, n);

      /* output */
      SCIPinfoMessage(scip, NULL, "\nFinal order of linear ordering constraint <%s>:\n", SCIPconsGetName(conss[c]));
      for (i = 0; i < n; ++i)
	 SCIPinfoMessage(scip, NULL, "%d ", indices[i]);
      SCIPinfoMessage(scip, NULL, "\n");

      SCIPfreeBufferArray(scip, &indices);
      SCIPfreeBufferArray(scip, &outdeg);
   }

   return SCIP_OKAY;
}

/** transforms constraint data into data belonging to the transformed problem */
static
SCIP_DECL_CONSTRANS(consTransLOP)
{  /*lint --e{715}*/
   SCIP_CONSDATA* consdata;
   SCIP_CONSDATA* sourcedata;
   int i;
   int j;
   int n;
   char s[SCIP_MAXSTRLEN];

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( sourcecons != NULL );
   assert( targetcons != NULL );

   SCIPdebugMsg(scip, "transforming linear ordering constraint <%s>.\n", SCIPconsGetName(sourcecons) );

   /* get data of original constraint */
   sourcedata = SCIPconsGetData(sourcecons);
   assert( sourcedata != NULL);

   /* create constraint data */
   SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );

   n = sourcedata->n;
   consdata->n = n;

   /* transform variables */
   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, n) );
   for (i = 0; i < n; ++i)
   {
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(consdata->vars[i]), n) ); /*lint !e866*/
      for (j = 0; j < n; ++j)
      {
	 if (j != i)
	 {
	    assert( sourcedata->vars[i][j] != NULL );
	    SCIP_CALL( SCIPgetTransformedVar(scip, sourcedata->vars[i][j], &(consdata->vars[i][j])) );
	 }
      }
   }

   /* create constraint */
   (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "t_%s", SCIPconsGetName(sourcecons));

   SCIP_CALL( SCIPcreateCons(scip, targetcons, s, conshdlr, consdata,
         SCIPconsIsInitial(sourcecons), SCIPconsIsSeparated(sourcecons),
         SCIPconsIsEnforced(sourcecons), SCIPconsIsChecked(sourcecons),
         SCIPconsIsPropagated(sourcecons), SCIPconsIsLocal(sourcecons),
         SCIPconsIsModifiable(sourcecons), SCIPconsIsDynamic(sourcecons),
         SCIPconsIsRemovable(sourcecons), SCIPconsIsStickingAtNode(sourcecons)) );

   return SCIP_OKAY;
}

/** LP initialization method of constraint handler */
static
SCIP_DECL_CONSINITLP(consInitlpLOP)
{  /*lint --e{715}*/
   char s[SCIP_MAXSTRLEN];
   int c;
   int nGen = 0;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( infeasible != NULL );

   *infeasible = FALSE;

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_VAR*** vars;
      int i;
      int j;
      int n;

      assert( conss != NULL );
      assert( conss[c] != NULL );
      SCIPdebugMsg(scip, "adding initial rows for linear ordering constraint <%s>.\n", SCIPconsGetName(conss[c]));

      consdata = SCIPconsGetData(conss[c]);
      assert( consdata != NULL );
      assert( consdata->vars != NULL );
      n = consdata->n;
      vars = consdata->vars;

      /* add symmetry equation */
      for (i = 0; i < n; ++i)
      {
	 for (j = i+1; j < n; ++j)
	 {
	    SCIP_ROW* row;

	    (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "sym#%d#%d", i, j);
	    SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, s, 1.0, 1.0, FALSE, FALSE, FALSE) );
	    SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
	    SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i][j], 1.0) );
	    SCIP_CALL( SCIPaddVarToRow(scip, row, vars[j][i], 1.0) );
	    SCIP_CALL( SCIPflushRowExtensions(scip, row) );
#ifdef SCIP_DEBUG
	    SCIPdebug( SCIPprintRow(scip, row, NULL) );
#endif
	    SCIP_CALL( SCIPaddRow(scip, row, FALSE, infeasible) );
	    SCIP_CALL( SCIPreleaseRow(scip, &row));
	    ++nGen;

            /* cannot handle infeasible case here - just exit */
            if ( *infeasible )
               return SCIP_OKAY;
	 }
      }
   }
   SCIPdebugMsg(scip, "added %d equations.\n", nGen);

   return SCIP_OKAY;
}

/** separation method of constraint handler for LP solutions */
static
SCIP_DECL_CONSSEPALP(consSepalpLOP)
{  /*lint --e{715}*/
   int nGen = 0;
   int c;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( conss != NULL );
   assert( result != NULL );

   *result = SCIP_DIDNOTRUN;

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_CONS* cons;
      SCIP_Bool cutoff;

      cons = conss[c];
      assert( cons != NULL );
      SCIPdebugMsg(scip, "separating LP solution for linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

      consdata = SCIPconsGetData(cons);
      assert( consdata != NULL );

      *result = SCIP_DIDNOTFIND;
      SCIP_CALL( LOPseparate(scip, conshdlr, consdata->n, consdata->vars, NULL, &nGen, &cutoff) );
      if ( cutoff )
      {
         *result = SCIP_CUTOFF;
         return SCIP_OKAY;
      }
   }
   if (nGen > 0)
      *result = SCIP_SEPARATED;
   SCIPdebugMsg(scip, "separated %d cuts.\n", nGen);

   return SCIP_OKAY;
}

/** separation method of constraint handler for arbitrary primal solutions */
static
SCIP_DECL_CONSSEPASOL(consSepasolLOP)
{  /*lint --e{715}*/
   int nGen = 0;
   int c;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( conss != NULL );
   assert( result != NULL );

   *result = SCIP_DIDNOTRUN;

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_CONS* cons;
      SCIP_Bool cutoff;

      cons = conss[c];
      assert( cons != NULL );
      SCIPdebugMsg(scip, "separating solution for linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

      consdata = SCIPconsGetData(cons);
      assert( consdata != NULL );

      *result = SCIP_DIDNOTFIND;
      SCIP_CALL( LOPseparate(scip, conshdlr, consdata->n, consdata->vars, sol, &nGen, &cutoff) );
      if ( cutoff )
      {
         *result = SCIP_CUTOFF;
         return SCIP_OKAY;
      }
   }
   if (nGen > 0)
      *result = SCIP_SEPARATED;

   return SCIP_OKAY;
}

/** constraint enforcing method of constraint handler for LP solutions */
static
SCIP_DECL_CONSENFOLP(consEnfolpLOP)
{  /*lint --e{715}*/
   char s[SCIP_MAXSTRLEN];
   int nGen = 0;
   int c;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( conss != NULL );
   assert( result != NULL );

   *result = SCIP_DIDNOTRUN;

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_CONS* cons;
      SCIP_VAR*** vars;
      int i;
      int j;
      int k;
      int n;

      cons = conss[c];
      assert( cons != NULL );
      SCIPdebugMsg(scip, "enforcing lp solution for linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

      consdata = SCIPconsGetData(cons);
      assert( consdata != NULL );

      n = consdata->n;
      vars = consdata->vars;
      assert( vars != NULL );

      for (i = 0; i < n; ++i)
      {
	 for (j = i + 1; j < n; ++j)
	 {
	    SCIP_Real valIJ;

	    valIJ = SCIPgetSolVal(scip, NULL, vars[i][j]);

	    /* if symmetry equations are violated - should not be the case, if they are added in the beginning */
	    if ( ! SCIPisFeasEQ(scip, 1.0 - valIJ, SCIPgetSolVal(scip, NULL, vars[j][i])) )
	    {
	       SCIP_ROW *row;
               SCIP_Bool infeasible;

	       (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "sym#%d#%d", i, j);

	       SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, s, 1.0, 1.0, FALSE, FALSE, TRUE) );
	       SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
	       SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i][j], 1.0) );
	       SCIP_CALL( SCIPaddVarToRow(scip, row, vars[j][i], 1.0) );
	       SCIP_CALL( SCIPflushRowExtensions(scip, row) );
#ifdef SCIP_DEBUG
	       SCIPdebug( SCIPprintRow(scip, row, NULL) );
#endif
	       SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
	       SCIP_CALL( SCIPreleaseRow(scip, &row));
	       ++nGen;

               if ( infeasible )
               {
                  *result = SCIP_CUTOFF;
                  return SCIP_OKAY;
               }
	    }

	    /* enforce triangle inequalities */
	    for (k = i + 1; k < n; ++k)
	    {
	       SCIP_Real sum;

	       if ( k == j )
		  continue;

	       sum = valIJ + SCIPgetSolVal(scip, NULL, vars[j][k]) + SCIPgetSolVal(scip, NULL, vars[k][i]);

	       /* if sum > 2.0, i.e., the cut is violated */
	       if ( SCIPisFeasGT(scip, sum, 2.0) ) /* this is the only difference to the separation call */
	       {
		  SCIP_ROW *row;
                  SCIP_Bool infeasible;

		  (void) SCIPsnprintf(s, SCIP_MAXSTRLEN, "triangle#%d#%d#%d", i, j, k);

		  SCIP_CALL( SCIPcreateEmptyRowConshdlr(scip, &row, conshdlr, s, -SCIPinfinity(scip), 2.0, FALSE, FALSE, TRUE) );
		  SCIP_CALL( SCIPcacheRowExtensions(scip, row) );
		  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[i][j], 1.0) );
		  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[j][k], 1.0) );
		  SCIP_CALL( SCIPaddVarToRow(scip, row, vars[k][i], 1.0) );
		  SCIP_CALL( SCIPflushRowExtensions(scip, row) );
#ifdef SCIP_DEBUG
		  SCIPdebug( SCIPprintRow(scip, row, NULL) );
#endif
		  SCIP_CALL( SCIPaddRow(scip, row, FALSE, &infeasible) );
		  SCIP_CALL( SCIPreleaseRow(scip, &row));
		  ++nGen;

                  if ( infeasible )
                  {
                     *result = SCIP_CUTOFF;
                     return SCIP_OKAY;
                  }
	       }
	    }
	 }
      }

      if (nGen > 0)
      {
	 *result = SCIP_SEPARATED;
	 return SCIP_OKAY;
      }

   }
   SCIPdebugMsg(scip, "all linear ordering constraints are feasible.\n");
   *result = SCIP_FEASIBLE;

   return SCIP_OKAY;
}

/** constraint enforcing method of constraint handler for pseudo solutions */
static
SCIP_DECL_CONSENFOPS(consEnfopsLOP)
{  /*lint --e{715}*/
   int c;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( conss != NULL );
   assert( result != NULL );

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_CONS* cons;
      SCIP_VAR*** vars;
      int i;
      int j;
      int k;
      int n;

      cons = conss[c];
      assert( cons != NULL );
      SCIPdebugMsg(scip, "enforcing pseudo solution for linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

      consdata = SCIPconsGetData(cons);
      assert( consdata != NULL );
      assert( consdata->vars != NULL );
      vars = consdata->vars;
      n = consdata->n;

      /* check triangle inequalities */
      for (i = 0; i < n; ++i)
      {
	 for (j = i + 1; j < n; ++j)
	 {
	    SCIP_Bool oneIJ;
            SCIP_Bool oneJI;

	    /* the priorities should ensure that the solution is integral */
	    assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, NULL, vars[i][j])) );
	    assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, NULL, vars[j][i])) );

            oneIJ = SCIPgetSolVal(scip, NULL, vars[i][j]) > 0.5 ? TRUE : FALSE;
            oneJI = SCIPgetSolVal(scip, NULL, vars[j][i]) > 0.5 ? TRUE : FALSE;

	    if ( oneIJ == oneJI )
	    {
	       SCIPdebugMsg(scip, "constraint <%s> infeasible (violated equation).\n", SCIPconsGetName(cons));
	       *result = SCIP_INFEASIBLE;
	       return SCIP_OKAY;
	    }

	    for (k = i + 1; k < n; ++k)
	    {
	       SCIP_Bool oneJK;
               SCIP_Bool oneKI;

	       if ( k == j )
		  continue;

	       assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, NULL, vars[j][k])) );
	       assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, NULL, vars[k][i])) );

	       oneJK = SCIPgetSolVal(scip, NULL, vars[j][k]) > 0.5 ? TRUE : FALSE;
	       oneKI = SCIPgetSolVal(scip, NULL, vars[k][i]) > 0.5 ? TRUE : FALSE;

	       /* if triangle inequality is violated */
	       if ( oneIJ && oneJK && oneKI )
	       {
		  SCIPdebugMsg(scip, "constraint <%s> infeasible (violated triangle ineq.).\n", SCIPconsGetName(cons));
		  *result = SCIP_INFEASIBLE;
		  return SCIP_OKAY;
	       }
	    }
	 }
      }
   }
   SCIPdebugMsg(scip, "all linear ordering constraints are feasible.\n");
   *result = SCIP_FEASIBLE;

   return SCIP_OKAY;
}

/** feasibility check method of constraint handler for integral solutions */
static
SCIP_DECL_CONSCHECK(consCheckLOP)
{  /*lint --e{715}*/
   int c;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( conss != NULL );
   assert( result != NULL );

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_CONS* cons;
      SCIP_VAR*** vars;
      int i;
      int j;
      int k;
      int n;

      cons = conss[c];
      assert( cons != NULL );
      SCIPdebugMsg(scip, "checking linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

      consdata = SCIPconsGetData(cons);
      assert( consdata != NULL );
      assert( consdata->vars != NULL );
      vars = consdata->vars;
      n = consdata->n;

      /* check triangle inequalities and symmetry equations */
      for (i = 0; i < n; ++i)
      {
	 for (j = i + 1; j < n; ++j)
	 {
	    SCIP_Bool oneIJ;
            SCIP_Bool oneJI;

	    /* the priorities should ensure that the solution is integral */
	    assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, sol, vars[i][j])) );
	    assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, sol, vars[j][i])) );

	    oneIJ = SCIPgetSolVal(scip, sol, vars[i][j]) > 0.5 ? TRUE : FALSE;
            oneJI = SCIPgetSolVal(scip, sol, vars[j][i]) > 0.5 ? TRUE : FALSE;

	    /* check symmetry equations */
	    if ( oneIJ == oneJI )
	    {
	       SCIPdebugMsg(scip, "constraint <%s> infeasible (violated equation).\n", SCIPconsGetName(cons));
	       *result = SCIP_INFEASIBLE;
               if( printreason )
               {
                  SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
                  SCIPinfoMessage(scip, NULL, "violation: symmetry equation violated <%s> = %.15g and <%s> = %.15g\n",
                     SCIPvarGetName(vars[i][j]), SCIPgetSolVal(scip, sol, vars[i][j]),
                     SCIPvarGetName(vars[j][i]), SCIPgetSolVal(scip, sol, vars[j][i]));
               }
	       return SCIP_OKAY;
	    }

	    for (k = i + 1; k < n; ++k)
	    {
	       SCIP_Bool oneJK;
               SCIP_Bool oneKI;

	       if ( k == j )
		  continue;

	       assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, sol, vars[j][k])) );
	       assert( SCIPisFeasIntegral(scip, SCIPgetSolVal(scip, sol, vars[k][i])) );

	       oneJK = SCIPisGT(scip, SCIPgetSolVal(scip, sol, vars[j][k]), 0.5);
	       oneKI = SCIPisGT(scip, SCIPgetSolVal(scip, sol, vars[k][i]), 0.5);

	       /* if triangle inequality is violated */
	       if ( oneIJ && oneJK && oneKI )
	       {
		  SCIPdebugMsg(scip, "constraint <%s> infeasible (violated triangle ineq.).\n", SCIPconsGetName(cons));
		  *result = SCIP_INFEASIBLE;
                  if( printreason )
                  {
                     SCIP_CALL( SCIPprintCons(scip, cons, NULL) );
                     SCIPinfoMessage(scip, NULL,
                        "violation: triangle inequality violated <%s> = %.15g, <%s> = %.15g, <%s> = %.15g\n",
                        SCIPvarGetName(vars[i][j]), SCIPgetSolVal(scip, sol, vars[i][j]),
                        SCIPvarGetName(vars[j][k]), SCIPgetSolVal(scip, sol, vars[j][k]),
                        SCIPvarGetName(vars[k][i]), SCIPgetSolVal(scip, sol, vars[k][i]));
                  }
		  return SCIP_OKAY;
	       }
	    }
	 }
      }
   }
   SCIPdebugMsg(scip, "all linear ordering constraints are feasible.\n");
   *result = SCIP_FEASIBLE;

   return SCIP_OKAY;
}

/** domain propagation method of constraint handler */
static
SCIP_DECL_CONSPROP(consPropLOP)
{  /*lint --e{715}*/
   int c;
   int nGen = 0;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( conss != NULL );
   assert( result != NULL );

   *result = SCIP_DIDNOTRUN;

   /* loop through all constraints */
   for (c = 0; c < nconss; ++c)
   {
      SCIP_CONSDATA* consdata;
      SCIP_CONS* cons;
      SCIP_VAR*** vars;
      int i;
      int j;
      int k;
      int n;

      cons = conss[c];
      assert( cons != NULL );
      SCIPdebugMsg(scip, "propagating linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

      *result = SCIP_DIDNOTFIND;

      consdata = SCIPconsGetData(cons);
      assert( consdata != NULL );
      assert( consdata->vars != NULL );

      vars = consdata->vars;
      n = consdata->n;

      /* check triangle inequalities */
      for (i = 0; i < n; ++i)
      {
	 for (j = i + 1; j < n; ++j)
	 {
            SCIP_Bool infeasible;
            SCIP_Bool tightened;

            /* if x[i][j] == 1 then x[j][i] = 0 */
	    if ( SCIPvarGetLbLocal(vars[i][j]) > 0.5 )
	    {
	       SCIP_CALL( SCIPinferBinvarCons(scip, vars[j][i], FALSE, cons, i*n + j, &infeasible, &tightened) );
	       if ( infeasible )
	       {
		  SCIPdebugMsg(scip, " -> node infeasible.\n");
                  SCIP_CALL( SCIPinitConflictAnalysis(scip, SCIP_CONFTYPE_PROPAGATION, FALSE) );
                  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[i][j]) );
                  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[j][i]) );
                  SCIP_CALL( SCIPanalyzeConflictCons(scip, cons, NULL) );
		  *result = SCIP_CUTOFF;
		  return SCIP_OKAY;
	       }
	       if ( tightened )
		  ++nGen;
	    }

	    /* if x[i][j] == 0 then x[j][i] = 1 */
	    if ( SCIPvarGetUbLocal(vars[i][j]) < 0.5 )
	    {
	       SCIP_CALL( SCIPinferBinvarCons(scip, vars[j][i], TRUE, cons, i*n + j, &infeasible, &tightened) );
	       if ( infeasible )
	       {
		  SCIPdebugMsg(scip, " -> node infeasible.\n");
                  SCIP_CALL( SCIPinitConflictAnalysis(scip, SCIP_CONFTYPE_PROPAGATION, FALSE) );
                  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[i][j]) );
                  SCIP_CALL( SCIPaddConflictBinvar(scip, vars[j][i]) );
                  SCIP_CALL( SCIPanalyzeConflictCons(scip, cons, NULL) );
		  *result = SCIP_CUTOFF;
		  return SCIP_OKAY;
	       }
	       if ( tightened )
		  ++nGen;
	    }

	    for (k = i + 1; k < n; ++k)
	    {
	       if ( k == j )
		  continue;

	       /* if x[i][j] == 1 and x[j][k] == 1 then x[k][i] = 0 */
	       if ( SCIPvarGetLbLocal(vars[i][j]) > 0.5 && SCIPvarGetLbLocal(vars[j][k]) > 0.5 )
	       {
		  SCIP_CALL( SCIPinferBinvarCons(scip, vars[k][i], FALSE, cons, n*n + i*n*n + j*n + k, &infeasible, &tightened) );
		  if ( infeasible )
		  {
		     SCIPdebugMsg(scip, " -> node infeasible.\n");
                     SCIP_CALL( SCIPinitConflictAnalysis(scip, SCIP_CONFTYPE_PROPAGATION, FALSE) );
                     SCIP_CALL( SCIPaddConflictBinvar(scip, vars[i][j]) );
                     SCIP_CALL( SCIPaddConflictBinvar(scip, vars[j][k]) );
                     SCIP_CALL( SCIPaddConflictBinvar(scip, vars[k][i]) );
                     SCIP_CALL( SCIPanalyzeConflictCons(scip, cons, NULL) );
		     *result = SCIP_CUTOFF;
		     return SCIP_OKAY;
		  }
		  if ( tightened )
		     ++nGen;
	       }

	       /* all other implications occur with other indices i, j, k */
	    }
	 }
      }
   }
   SCIPdebugMsg(scip, "propagated %d domains.\n", nGen);
   if (nGen > 0)
      *result = SCIP_REDUCEDDOM;

   return SCIP_OKAY;
}

/** propagation conflict resolving method of constraint handler */
static
SCIP_DECL_CONSRESPROP(consRespropLOP)
{  /*lint --e{715}*/
   SCIP_CONSDATA* consdata;
   SCIP_VAR*** vars;
   int n;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( cons != NULL );
   assert( infervar != NULL );
   assert( bdchgidx != NULL );
   assert( result != NULL );

   SCIPdebugMsg(scip, "Propagation resolution of constraint <%s>.\n", SCIPconsGetName(cons));
   *result = SCIP_DIDNOTFIND;

   consdata = SCIPconsGetData(cons);
   assert( consdata != NULL);
   assert( consdata->vars != NULL );

   n = consdata->n;
   vars = consdata->vars;

   assert( 0 <= inferinfo && inferinfo < n*n + n*n*n );

   /* if the conflict came from an equation */
   if ( inferinfo < (n*n) )
   {
      int index1;
      int index2;

      index1 = inferinfo/n;
      index2 = inferinfo % n;
      assert( 0 <= index1 && index1 < n );
      assert( 0 <= index2 && index2 < n );
      assert( vars[index2][index1] == infervar );

      /* if the variable was fixed to 0 */
      if ( SCIPvarGetUbAtIndex(infervar, bdchgidx, FALSE) > 0.5 && SCIPvarGetUbAtIndex(infervar, bdchgidx, TRUE) < 0.5 )
      {
	 SCIPdebugMsg(scip, " -> reason for x[%d][%d] == 0 was x[%d][%d] = 1.\n", index2, index1, index1, index2);
	 /* the reason was that x[i][j] was fixed to 1 */
	 SCIP_CALL( SCIPaddConflictLb(scip, vars[index1][index2], bdchgidx) );
	 *result = SCIP_SUCCESS;
	 return SCIP_OKAY;
      }

      /* if the variable was fixed to 1 */
      if ( SCIPvarGetLbAtIndex(infervar, bdchgidx, FALSE) < 0.5 && SCIPvarGetLbAtIndex(infervar, bdchgidx, TRUE) > 0.5 )
      {
	 SCIPdebugMsg(scip, " -> reason for x[%d][%d] == 1 was x[%d][%d] = 0.\n", index2, index1, index1, index2);
	 /* the reason was that x[i][j] was fixed to 0 */
	 SCIP_CALL( SCIPaddConflictUb(scip, vars[index1][index2], bdchgidx) );
	 *result = SCIP_SUCCESS;
	 return SCIP_OKAY;
      }
   }
   else
   {
      /* otherwise the conflict came from a triangle inequality */
      int index1;
      int index2;
      int index3;

      index1 = (inferinfo - n*n)/(n*n);
      index2 = (inferinfo - n*n - index1 * n*n)/n;
      index3 = (inferinfo - n*n) % n;

      assert( 0 <= index1 && index1 < n );
      assert( 0 <= index2 && index2 < n );
      assert( 0 <= index3 && index3 < n );
      assert( index1 != index2 && index2 != index3 && index1 != index3 );
      assert( vars[index3][index1] == infervar );

      /* the variable should have been fixed to 0 */
      assert( SCIPvarGetUbAtIndex(infervar, bdchgidx, FALSE) > 0.5 && SCIPvarGetUbAtIndex(infervar, bdchgidx, TRUE) < 0.5 );

      /* the reason was that x[index1][index2] and x[index2][index3] were fixed to 1 */
      SCIPdebugMsg(scip, " -> reason for x[%d][%d] == 0 was x[%d][%d] = x[%d][%d] = 1.\n", index3, index1, index1, index2, index2, index3);
      SCIP_CALL( SCIPaddConflictLb(scip, vars[index1][index2], bdchgidx) );
      SCIP_CALL( SCIPaddConflictLb(scip, vars[index2][index3], bdchgidx) );
      *result = SCIP_SUCCESS;
   }

   return SCIP_OKAY;
}

/** variable rounding lock method of constraint handler */
static
SCIP_DECL_CONSLOCK(consLockLOP)
{  /*lint --e{715}*/
   int i;
   int j;
   SCIP_CONSDATA* consdata;
   SCIP_VAR*** vars;
   int n;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( cons != NULL );

   SCIPdebugMsg(scip, "Locking linear ordering constraint <%s>.\n", SCIPconsGetName(cons));

   /* get data of constraint */
   consdata = SCIPconsGetData(cons);
   assert( consdata != NULL);
   assert( consdata->vars != NULL );
   n = consdata->n;
   vars = consdata->vars;

   for (i = 0; i < n; ++i)
   {
      for (j = 0; j < n; ++j)
      {
	 if ( i != j )
	 {
	    /* the constraint may be violated in any way */
	    SCIP_CALL( SCIPaddVarLocksType(scip, vars[i][j], SCIP_LOCKTYPE_MODEL, nlockspos + nlocksneg, nlockspos + nlocksneg) );
	 }
      }
   }

   return SCIP_OKAY;
}

/** constraint display method of constraint handler */
static
SCIP_DECL_CONSPRINT(consPrintLOP)
{  /*lint --e{715}*/
   SCIP_CONSDATA* consdata;
   SCIP_VAR*** vars;
   int i;
   int j;
   int n;

   assert( scip != NULL );
   assert( conshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(conshdlr), CONSHDLR_NAME) == 0 );
   assert( cons != NULL );

   consdata = SCIPconsGetData(cons);
   assert( consdata != NULL );
   assert( consdata->vars != NULL );
   n = consdata->n;
   vars = consdata->vars;

   SCIPinfoMessage(scip, file, "LOP[");
   for (i = 0; i < n; ++i)
   {
      if ( i > 0 )
	 SCIPinfoMessage(scip, file, ", ");
      SCIPinfoMessage(scip, file, "(");
      for (j = 0; j < n; ++j)
      {
	 if ( j != i )
	 {
	    if ( j > 0 && (i > 0 || j > 1) )
	       SCIPinfoMessage(scip, file, ",");
	    SCIPinfoMessage(scip, file, "%s", SCIPvarGetName(vars[i][j]));
	 }
      }
      SCIPinfoMessage(scip, file, ")");
   }
   SCIPinfoMessage(scip, file, "]\n");

   return SCIP_OKAY;
}

/** constraint copying method of constraint handler */
static
SCIP_DECL_CONSCOPY(consCopyLOP)
{  /*lint --e{715}*/
   SCIP_CONSDATA* sourcedata;
   SCIP_VAR*** sourcevars;
   SCIP_VAR*** vars;
   int i;
   int j;
   int n;

   assert( scip != NULL );
   assert( sourceconshdlr != NULL );
   assert( strcmp(SCIPconshdlrGetName(sourceconshdlr), CONSHDLR_NAME) == 0 );
   assert( cons != NULL );
   assert( sourcescip != NULL );
   assert( sourcecons != NULL );
   assert( varmap != NULL );
   assert( valid != NULL );

   *valid = TRUE;

   SCIPdebugMsg(scip, "Copying method for linear ordering constraint handler.\n");

   sourcedata = SCIPconsGetData(sourcecons);
   assert( sourcedata != NULL );

   n = sourcedata->n;
   sourcevars = sourcedata->vars;
   assert( sourcevars != NULL );

   SCIP_CALL( SCIPallocBufferArray(scip, &vars, n) );
   BMSclearMemoryArray(vars, n);

   for (i = 0; i < n; ++i)
   {
      SCIP_CALL( SCIPallocBufferArray(scip, &(vars[i]), n) ); /*lint !e866*/

      for (j = 0; j < n && *valid; ++j)
      {
         if ( i != j )
         {
            SCIP_CALL( SCIPgetVarCopy(sourcescip, scip, sourcevars[i][j], &vars[i][j], varmap, consmap, global, valid) );
            assert( !(*valid) || vars[i][j] != NULL );
         }
      }
   }

   if ( *valid )
   {
      /* create copied constraint */
      if ( name == 0 )
         name = SCIPconsGetName(sourcecons);

      SCIP_CALL( SCIPcreateConsLOP(scip, cons, name, n, vars,
            initial, separate, enforce, check, propagate, local, modifiable, dynamic, removable, stickingatnode) );
   }

   /* free memory in reverse order */
   for (i = n-1; i >= 0; --i)
      SCIPfreeBufferArrayNull(scip, &vars[i]);
   SCIPfreeBufferArray(scip, &vars);

   return SCIP_OKAY;
}

/** creates the handler for linear ordering constraints and includes it in SCIP */
SCIP_RETCODE SCIPincludeConshdlrLOP(
   SCIP*                 scip                /**< SCIP data structure */
   )
{
   SCIP_CONSHDLR* conshdlr = NULL;

   /* include constraint handler */
   SCIP_CALL( SCIPincludeConshdlrBasic(scip, &conshdlr, CONSHDLR_NAME, CONSHDLR_DESC,
         CONSHDLR_ENFOPRIORITY, CONSHDLR_CHECKPRIORITY, CONSHDLR_EAGERFREQ, CONSHDLR_NEEDSCONS,
         consEnfolpLOP, consEnfopsLOP, consCheckLOP, consLockLOP, NULL) );
   assert( conshdlr != NULL );

   SCIP_CALL( SCIPsetConshdlrDelete(scip, conshdlr, consDeleteLOP) );
   SCIP_CALL( SCIPsetConshdlrExit(scip, conshdlr, consExitLOP) );
   SCIP_CALL( SCIPsetConshdlrCopy(scip, conshdlr, conshdlrCopyLOP, consCopyLOP) );
   SCIP_CALL( SCIPsetConshdlrTrans(scip, conshdlr, consTransLOP) );
   SCIP_CALL( SCIPsetConshdlrInitlp(scip, conshdlr, consInitlpLOP) );
   SCIP_CALL( SCIPsetConshdlrSepa(scip, conshdlr, consSepalpLOP, consSepasolLOP,
         CONSHDLR_SEPAFREQ, CONSHDLR_SEPAPRIORITY, CONSHDLR_DELAYSEPA) );
   SCIP_CALL( SCIPsetConshdlrProp(scip, conshdlr, consPropLOP, CONSHDLR_PROPFREQ,
         CONSHDLR_DELAYPROP, CONSHDLR_PROP_TIMING) );
   SCIP_CALL( SCIPsetConshdlrResprop(scip, conshdlr, consRespropLOP) );
   SCIP_CALL( SCIPsetConshdlrPrint(scip, conshdlr, consPrintLOP) );

   return SCIP_OKAY;
}

/** creates and captures a linear ordering constraint */
SCIP_RETCODE SCIPcreateConsLOP(
   SCIP*                 scip,               /**< SCIP data structure */
   SCIP_CONS**           cons,               /**< pointer to hold the created constraint */
   const char*           name,               /**< name of constraint */
   int                   n,                  /**< number of elements */
   SCIP_VAR***           vars,               /**< n x n matrix of binary variables */
   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 subject to aging? */
   SCIP_Bool             removable,          /**< should the relaxation be removed from the LP due to aging or cleanup? */
   SCIP_Bool             stickingatnode      /**< should the constraint always be kept at the node where it was added, even
                                              *   if it may be moved to a more global node? */
   )
{
   SCIP_CONSHDLR* conshdlr;
   SCIP_CONSDATA* consdata;
   int i;
   int j;

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

   /* create constraint data */
   SCIP_CALL( SCIPallocBlockMemory(scip, &consdata) );

   consdata->n = n;
   SCIP_CALL( SCIPallocBlockMemoryArray(scip, &consdata->vars, n) );
   for (i = 0; i < n; ++i)
   {
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &(consdata->vars[i]), n) ); /*lint !e866*/
      for (j = 0; j < n; ++j)
      {
	 if ( j != i )
	 {
	    assert( vars[i][j] != NULL );
	    consdata->vars[i][j] = vars[i][j];
	 }
      }
   }

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

   return SCIP_OKAY;
}
