
#include "handle_localcuts.h"


#define MAXNSOLS               10    /**< maximal number of additional solutions allowed for oracle call */
#ifdef SCIP_DEBUG
#define OUTPUTFREQ              1    /**< frequency for additional output */
#endif
#define MIN_SCALEVAL       1.0e-4    /**< minimal scaling value allowed for cuts to be generated */
#define MAXINTEGER        10000.0    /**< maximal absolute value of scaling factor for making cuts integral */
#define MAXDENOMINATOR    10000LL    /**< maximal value of the denominator in rational approximation of cuts */
#define MAXABSCOEF        10000LL    /**< maximal absolute value of a coefficient in a local cut to be accepted */

#define INIT_NMAXPOINTS      1000    /**< initial maximal number of points */


/* ------------ Basic local cut methods ---------------------------------------------------------- */

/** Finds rational approximation to given real number.
 *  @return value of (signed) error
 *
 * (original version by David Eppstein / UC Irvine / 8 Aug 1993)
 * (slight modifications by Giovanni Rinaldi / August 2010)
 *
 *  Based on the theory of continued fractions. If
 *  \f[
 *     x = a_1 + \frac{1}{a_2 + \frac{1}{a_3 + \frac{1}{a_4 + \dots}}},
 *  \f]
 *  then best approximation is found by truncating this series (with some adjustments in the last term).
 *
 *  Note the fraction can be recovered as the first column of the matrix
 *  \f[
 *     \left(\begin{array}{rr}
 *       a_1 & 1\\
 *         1 & 0
 *     \end{array}\right)
 *     \cdot
 *     \left(\begin{array}{rr}
 *       a_2 & 1\\
 *         1 & 0
 *     \end{array}\right)
 *     \cdot
 *     \left(\begin{array}{rr}
 *       a_3 & 1\\
 *         1 & 0
 *     \end{array}\right)
 *     \cdot \dots
 *  \f]
 *
 *  Instead of keeping the sequence of continued fraction terms, we just keep the last partial product of these
 *  matrices.
 */
static
SCIP_Real rationalApproximation(
   SCIP_Real             x,                  /**< real number to approximate */
   SCIP_Longint          maxden,             /**< maximum allowed denominator */
   SCIP_Longint*         num,                /**< returned numerator */
   SCIP_Longint*         den,                /**< returned denominator */
   int                   appx                /**< desired approximation direction (1: upper, -1: lower, 0: best);
                                              *   if the error is below @p tol this parameter has no effect. */
   )
{
   SCIP_Real startx;
   SCIP_Longint n1;
   SCIP_Longint d1;
   SCIP_Real e1;
   SCIP_Longint n2;
   SCIP_Longint d2;
   SCIP_Real e2;
   SCIP_Real tol = (SCIP_Real) (1.0 / SCIP_LONGINT_MAX);
   int c1;
   int sign;
   SCIP_Longint m[2][2];
   SCIP_Longint ai;

   sign = x > 0.0 ? 1 : -1;
   startx = REALABS(x);

   if ( startx > SCIP_LONGINT_MAX )
   {
      SCIPerrorMessage("Input number too large!\n");
      goto failure;
   }

   m[0][0] = m[1][1] = 1;
   m[0][1] = m[1][0] = 0;
   ai = (SCIP_Longint) x;
   while ( ABS(m[1][0] * ai + m[1][1]) <= maxden )
   {
      SCIP_Real err;
      SCIP_Longint t;

      t = m[0][0] * ai + m[0][1];
      m[0][1] = m[0][0];
      m[0][0] = t;

      t = m[1][0] * ai + m[1][1];
      m[1][1] = m[1][0];
      m[1][0] = t;

      x -= (SCIP_Real) ai;
      err = REALABS((SCIP_Real) m[0][0] / (SCIP_Real) m[1][0]) - startx;

      if ( REALABS(err) < tol )
         break;

      if ( REALABS(x) < tol )
         break;

      x = 1.0 / x;

      ai = (SCIP_Longint) x;
   }

   n1 = m[0][0];
   d1 = m[1][0];
   e1 = (SCIP_Real) sign * (REALABS((SCIP_Real) n1 / (SCIP_Real) d1) - startx);

   ai = (maxden - m[1][1]) / m[1][0]; /* division in integers */
   n2 = m[0][0] * ai + m[0][1];
   d2 = m[1][0] * ai + m[1][1];
   e2 = (SCIP_Real) sign * (REALABS((SCIP_Real) n2 / (SCIP_Real) d2) - startx);

   if ( appx > 0 )
      c1 = (e1 > -tol);
   else if ( appx < 0 )
      c1 = (e1 < tol);
   else
      c1 = REALABS(e1) < REALABS(e2);

   if ( c1 )
   {
      *num = sign * ABS(n1);
      *den = ABS(d1);
      return (e1);
   }
   else
   {
      *num = sign * ABS(n2);
      *den = ABS(d2);
      return (e2);
   }

 failure:
   *num = (SCIP_Longint) (sign * trunc(startx));
   *den = 1;

   return (*num - sign * startx);
}


/** Round coefficient vector to integer values */
static
SCIP_RETCODE makeInteger(
   SCIP*                 scip,               /**< SCIP pointer */
   int                   dim,                /**< dimension of vector */
   const SCIP_Real*      x,                  /**< vector */
   SCIP_Real             maxint,             /**< maximal integer value allowed for scaling */
   SCIP_Longint          maxden,             /**< maximal denominator allowed */
   SCIP_Real*            scale               /**< scaling value to make x integral */
   )
{
   int i;
   int s;
   SCIP_Longint g;
   int udim;
   SCIP_Real* u;
   SCIP_Real lastu;
   SCIP_Longint p;
   SCIP_Longint num;
   SCIP_Longint den;

   assert( x != NULL );
   assert( scale != NULL );

   /* We assume that 1/maxint is above the feasibility tolerance, i.e., maxint <= 1e6 for default settings; otherwise we
    * need to turn the assert for zero coefficients below into a check. */
   assert( 1.0/maxint >= SCIPfeastol(scip) );

   /* store absolute values of x in u */
   SCIP_CALL( SCIPallocBufferArray(scip, &u, dim) );
   for (i = 0, udim = 0; i < dim; ++i)
   {
      SCIP_Real val;

      val = REALABS(x[i]);

      /* take entry only if it is large enough */
      if ( val * maxint + 0.5 >= 1.0 )
      {
         assert( ! SCIPisFeasZero(scip, val) );
         u[udim++] = val;
      }
   }

   if ( udim == 0 )
   {
      *scale = 0;
      SCIPfreeBufferArray(scip, &u);
      return SCIP_OKAY;
   }

   /* sort vector u non-decreasingly */
   SCIPsortReal(u, udim);

   /* extract s distinct strictly increasing values in u */
   for (i = 1, s = 0; i < udim; ++i)
   {
      if ( SCIPisFeasGT(scip, u[i], u[s]) )
         u[++s] = u[i];
   }
   ++s;

   /* divide by the first number */
   for (i = 1; i < s; ++i)
      u[i] /= u[0];
   *scale = 1.0/u[0];
   u[0] = 1.0;

   /* initial scaling value */
   p = 1LL;

   /* loop through increasing values */
   lastu = 1.0;
   for (i = 1; i < s; ++i)
   {
      /* if value is too large */
      if ( u[i] * (SCIP_Real) p > maxint )
      {
         *scale = 0;
         SCIPfreeBufferArray(scip, &u);
         return SCIP_OKAY;
      }

      den = 1;
      lastu = u[i];
      if (! SCIPisIntegral(scip, u[i] * (SCIP_Real) p) )
      {
         (void) rationalApproximation(u[i], maxden, &num, &den, 0);
#if 0
         u[i] = (SCIP_Real) num; /* use later */
#endif

         /* g = getGcd (p, den); */
         g = SCIPcalcGreComDiv(p, den);
         assert( den % g == 0 );
         den /= g;
         p *= den;
      }
   }
   assert( s >= 1 );

   if ( lastu * (SCIP_Real) p > maxint )
      *scale = 0.0;
   else
      *scale *= (SCIP_Real) p;

   SCIPfreeBufferArray(scip, &u);

   return SCIP_OKAY;
}


/** check and create cut */
static
SCIP_RETCODE LC_checkCut(
   SCIP*                 scip,               /**< SCIP pointer */
   SCIP_Bool             nicecut,            /**< try to make a nice cut by finding a rational approximation */
   const SCIP_Real*      x,                  /**< point to be separated in reduced space */
   int                   dim,                /**< dimension of subproblem */
   const SCIP_Real*      coeff,              /**< coefficients of the inequality */
   SCIP_Real             maxint,             /**< maximal integer value allowed for scaling */
   SCIP_Longint          maxden,             /**< maximal denominator allowed */
   SCIP_Longint          maxabscoef,         /**< maximal absolute value of a coefficient in a cut to be accepted */
   LC_ORACLE((*Oracle)),                     /**< optimization oracle callback */
   const LC_DATA*        data,               /**< data passed to oracle */
   SCIP_Real**           oraclesols,         /**< temporary vector for oracle solution */
   SCIP_Real*            oracleobjs,         /**< temporary vector for oracle objectives */
   SCIP_Real*            cutcoefs,           /**< cut coefficients of (tilted) cut (or NULL) */
   SCIP_Real*            cutrhs,             /**< pointer to store the rhs of the cut (or NULL) */
   SCIP_Bool*            success             /**< pointer to store whether cut generation was successful */
   )
{
   LC_STATUS status;
   SCIP_Real sidevalue;
   SCIP_Real rhs;
   SCIP_Real sum = 0.0;
   int noraclesols;
   int cnt = 0;
   int j;

   *success = FALSE;

   /* compute scaling to get integer vector */
   if ( nicecut )
   {
      SCIP_Real scale = 0.0;
      int lpdim;

      lpdim = dim + 1;

      SCIP_CALL( makeInteger(scip, lpdim, coeff, maxint, maxden, &scale) );

      if ( MIN_SCALEVAL < scale )
      {
         /* the following are all reals because this has to be passed to the oracle */
         SCIP_Real* Coeff;
         SCIP_Longint g = 1;
         SCIP_Longint maxabs = -1;

         /* vector of reduced coefficients - dense because of oracle */
         SCIP_CALL( SCIPallocBufferArray(scip, &Coeff, lpdim) );

         /* round coefficient vector and compute gcd */
         for (j = 0; j < dim; ++j)
         {
            if ( ! SCIPisFeasZero(scip, coeff[j]) )
            {
               assert( SCIPisFeasGE(scip, REALABS(coeff[j] * scale), 1.0) );
               Coeff[j] = SCIPround(scip, coeff[j] * scale);
               assert( SCIPisIntegral(scip, Coeff[j]) );
               if ( cnt == 0 )
                  g = (SCIP_Longint) REALABS(Coeff[j]);
               else
               {
                  if ( g != 1 )
                     g = SCIPcalcGreComDiv(g, (SCIP_Longint) REALABS(Coeff[j]));
               }
               ++cnt;
            }
            else
               Coeff[j] = 0.0;
         }
         assert( g > 0 );

         rhs = SCIPround(scip, coeff[dim] * scale);
         Coeff[dim] = rhs;

         /* divide by gcd and compute sum */
         cnt = 0;
         for (j = 0; j < dim; ++j)
         {
            if ( ! SCIPisZero(scip, Coeff[j]) )
            {
               if ( g != 1 )
                  Coeff[j] /= (SCIP_Real) g;

               if ( REALABS(Coeff[j]) > maxabs )
                  maxabs = (SCIP_Longint) REALABS(Coeff[j]);

               sum += x[j] * Coeff[j];
               cutcoefs[j] = Coeff[j];
            }
            else
               cutcoefs[j] = 0.0;
         }

         if ( g != 1 )
            rhs /= (SCIP_Real) g;

         /* only consider cuts with absolute values of coefficients at most maxabscoef */
         if ( maxabs <= maxabscoef )
         {
            /* call oracle with rounded coefficent vector as objective */
            SCIP_CALL( Oracle(scip, dim, data, Coeff, -SCIPinfinity(scip), 1, &noraclesols, oraclesols, oracleobjs, &status) );

            *cutrhs = oracleobjs[0];

            /* only proceed if we solved the problem */
            if ( status == LC_STATUS_OPTIMAL )
            {
               /* check whether optimal solution value remains the same (for target cuts just compute the value) */
               sidevalue = oracleobjs[0];
               if ( ! SCIPisEQ(scip, sidevalue, rhs) )
               {
                  SCIPwarningMessage(scip, "Rhs had to be changed! (new: %f, computed: %f)\n", sidevalue, rhs);
               }

               /* create cut if it is violated and valid in the case of target cuts */
               if ( SCIPisEfficacious(scip, sum - sidevalue) )
               {
                  *success = TRUE;
               }
               else
               {
                  SCIPdebugMsg(scip, "Cut is not efficious (violation: %f).\n", sidevalue - sum);
               }
            }
         }
         else
         {
            SCIPdebugMsg(scip, "Maximal absolute value of cut %" SCIP_LONGINT_FORMAT " larger than bound %" SCIP_LONGINT_FORMAT ".\n", maxabs, maxabscoef);
         }
         SCIPfreeBufferArray(scip, &Coeff);
      }
      else
      {
         SCIPdebugMsg(scip, "Coefficients get too large: inequality rejected!\n");
      }
   }
   else
   {
      /* call oracle with coefficent vector as objective */
      SCIP_CALL( Oracle(scip, dim, data, coeff, -SCIPinfinity(scip), 1, &noraclesols, oraclesols, oracleobjs, &status) );

      rhs = coeff[dim];

      /* only proceed if we solved the problem */
      if ( status == LC_STATUS_OPTIMAL )
      {
         /* check whether optimal solution value remains the same (for target cuts just compute the value) */
         sidevalue = oracleobjs[0];
         if ( ! SCIPisEQ(scip, sidevalue, rhs) )
         {
            SCIPwarningMessage(scip, "Rhs had to be changed! (new: %f, computed: %f)\n", sidevalue, rhs);
         }
         *cutrhs = sidevalue;

         assert( cnt == 0 );
         for (j = 0; j < dim; ++j)
         {
            if ( ! SCIPisZero(scip, coeff[j]) )
            {
               sum += x[j] * coeff[j];
               cutcoefs[j] = coeff[j];
            }
            else
               cutcoefs[j] = 0.0;
         }

         /* create cut if it is violated and valid in the case of target cuts */
         if ( SCIPisEfficacious(scip, sum - sidevalue) )
         {
            *success = TRUE;
         }
         else
         {
            SCIPdebugMsg(scip, "Cut is not efficious (violation: %f).\n", sidevalue - sum);
         }
      }
   }

   return SCIP_OKAY;
}


/** compute relative violation \f$\frac{a^T x - a_0}{a^T(x - x^0)}\f$ for cut \f$a^T x \leq a_0\f$, point \f$x\f$ and relative interior point \f$x^0\f$ */
static
SCIP_Real computeGamma(
   const int             dim,                /**< dimension of problem */
   const SCIP_Real*      x,                  /**< point to be separated */
   const SCIP_Real*      xi,                 /**< relative interior point */
   const SCIP_Real*      coeff               /**< cut coefficients */
   )
{
  SCIP_Real ce = 0.0;
  SCIP_Real cx = 0.0;
  SCIP_Real den;
  SCIP_Real num;
  int i;

  for (i = 0; i < dim; ++i)
  {
     ce += xi[i] * coeff[i];
     cx += x[i] * coeff[i];
  }

  den = cx - ce;
  num = cx - coeff[dim];

  if ( den < -1.0e-8 || den > 1.0e-8 )
     return (num / den);

  return (-1.0);
}

/** perturb cut */
static
SCIP_RETCODE perturb(
   SCIP*                 scip,               /**< SCIP pointer */
   SCIP_LPI*             masterlp,           /**< LP */
   int                   dim,                /**< dimension of problem */
   const SCIP_Real*      x,                  /**< current point to be separated */
   const SCIP_Real*      xi,                 /**< relative interior point */
   SCIP_Real*            coeff,              /**< coefficients of cut */
   SCIP_Bool*            perturbed           /**< pointer to store whether perturbation was performed */
   )
{
   SCIP_Real* rhs = NULL;
   SCIP_Real* lhs = NULL;
   SCIP_Real g;
   int* rind = NULL;
   int i;

   assert( scip != NULL );
   assert( masterlp != NULL );
   assert( x != NULL );
   assert( xi != NULL );
   assert( coeff != NULL );
   assert( perturbed != NULL );

   *perturbed = FALSE;

   /* compute relative violation */
   g = computeGamma(dim, x, xi, coeff);

   if ( g >= 0.0 )
   {
      SCIP_CALL( SCIPallocBufferArray(scip, &rind, dim) );
      SCIP_CALL( SCIPallocBufferArray(scip, &rhs, dim) );
      SCIP_CALL( SCIPallocBufferArray(scip, &lhs, dim) );

      for (i = 0; i < dim; ++i)
      {
         rind[i] = i;
         lhs[i] = -SCIPlpiInfinity(masterlp);
         rhs[i] = x[i] * (1.0 - g) + xi[i] * g;
      }

      SCIP_CALL( SCIPlpiChgSides(masterlp, dim, rind, lhs, rhs) );

      *perturbed = TRUE;

      SCIPfreeBufferArray(scip, &rind);
      SCIPfreeBufferArray(scip, &rhs);
      SCIPfreeBufferArray(scip, &lhs);
   }

   return SCIP_OKAY;
}


/** ensures space for points */
static
SCIP_RETCODE ensurePointNumber(
   SCIP*                 scip,               /**< SCIP data structure */
   int**                 pointbeg,           /**< beginning of points */
   int*                  num,                /**< current size */
   int                   newnum              /**< minimum number of entries to store */
   )
{
   assert( scip != NULL );
   assert( pointbeg != NULL );
   assert( num != NULL );

   if ( newnum > *num )
   {
      int newsize;

      newsize = SCIPcalcMemGrowSize(scip, newnum);
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, pointbeg, *num, newsize) );
      *num = newsize;
   }

   return SCIP_OKAY;
}


/** ensures space for points nonzeros */
static
SCIP_RETCODE ensurePointSize(
   SCIP*                 scip,               /**< SCIP data structure */
   int**                 pointind,           /**< nonzero component indices */
   SCIP_Real**           pointval,           /**< nonzero values */
   int*                  num,                /**< current size */
   int                   newnum              /**< minimum number of entries to store */
   )
{
   assert( scip != NULL );
   assert( pointind != NULL );
   assert( pointval != NULL );
   assert( num != NULL );

   if ( newnum > *num )
   {
      int newsize;

      newsize = SCIPcalcMemGrowSize(scip, newnum);
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, pointind, *num, newsize) );
      SCIP_CALL( SCIPreallocBlockMemoryArray(scip, pointval, *num, newsize) );
      *num = newsize;
   }

   return SCIP_OKAY;
}


//------------ local cut methods -----------------------------------------------------------------


#ifdef SCIP_CHECK
/** check whether point is new (for debugging only) */
static
void checkPointIsNew(
   SCIP*                 scip,               /**< SCIP pointer */
   SCIP_Real*            point,              /**< new point? */
   int                   dim,                /**< dimension of point */
   int                   npoints,            /**< current number of points */
   int*                  pointbeg,           /**< index of points in pointind/pointval arrays */
   int*                  pointind,           /**< indices of nonzero components of created points */
   SCIP_Real*            pointval            /**< values of nonzero components of created points */
)
{
   int i;
   int k;
   int j;
   int ind;

   /* loop through all points */
   for (i = 0; i < npoints; ++i)
   {
      SCIP_Bool equal = TRUE;

      /* check nonzeros of stored points */
      for (j = pointbeg[i]; j < pointbeg[i+1]; ++j)
      {
         ind = pointind[j];
         assert( 0 <= ind && ind < dim );
         if ( ! SCIPisFeasEQ(scip, pointval[j], point[ind]) )
         {
            equal = FALSE;
            break;
         }
      }

      /* check nonzeros of point */
      if ( equal )
      {
         k = 0;
         j = pointbeg[i];
         while ( k < dim && j < pointbeg[i+1] )
         {
            while ( k < dim && SCIPisFeasZero(scip, point[k]) )
               ++k;

            if ( k < dim )
            {
               /* use that pointind is sorted! */
               while ( j < pointbeg[i+1] && pointind[j] < k )
                  ++j;

               if ( j < pointbeg[i+1] )
               {
                  if ( pointind[j] > k )
                  {
                     equal = FALSE;
                     break;
                  }

                  assert( pointind[j] == k );
                  if ( ! SCIPisFeasEQ(scip, pointval[j], point[k]) )
                  {
                     equal = FALSE;
                     break;
                  }
                  ++k;
                  ++j;
               }
            }
         }
         /* find final nonzero - if necessary */
         while ( k < dim && SCIPisFeasZero(scip, point[k]) )
            ++k;
         if ( j >= pointbeg[i+1] && k < dim )
            equal = FALSE;
      }
      if ( equal )
      {
         printf("The new point:\n");
         for (k = 0; k < dim; ++k)
         {
            if ( ! SCIPisZero(scip, point[k]) )
               printf("%f [%d] ", point[k], k);
         }
         printf("\nis equal to the old point (%d | %d)\n", i, npoints);
         for (j = pointbeg[i]; j < pointbeg[i+1]; ++j)
            printf("%f [%d] ", pointval[j], pointind[j]);
         printf("\n");
      }
      assert( ! equal );
   }
}
#endif


/** check whether point is new (for debugging only) */
static
SCIP_RETCODE checkPoint(
   SCIP*                 scip,               /**< SCIP pointer */
   SCIP_Real*            point,              /**< new point? */
   int                   dim,                /**< dimension of point */
   SCIP_VAR**            vars                /**< original variables */
   )
{
   SCIP_Bool stored;
   SCIP_SOL* sol;
   int j;

   SCIP_CALL( SCIPcreateSol(scip, &sol, NULL) );

   for (j = 0; j < dim; ++j)
   {
      if ( ! SCIPisFeasZero(scip, point[j]) )
      {
         SCIP_CALL( SCIPsetSolVal(scip, sol, vars[j], point[j]) );
      }
   }

   SCIP_CALL( SCIPtrySolFree(scip, &sol, FALSE, FALSE, TRUE, TRUE, TRUE, &stored) );

#ifdef SCIP_DEBUG
   if ( stored )
      SCIPdebugMsg(scip, "Found solution!\n");
#endif

   return SCIP_OKAY;
}


/** Generates a single local cut
 *
 * - The procedure finds an local cut violated by the point stored in \f$x_1, \ldots,
 *   x_{\mbox{dim}}\f$.
 * - The cuts are added as rows to the current LP.
 * - It uses a time limit of @p maxtime seconds of CPU time (if @p maxtime is not
 *   NULL). The total used CPU time in seconds is returned in @p usedtime.
 * - The point @p xrelint should be in the relative interior (needed for tilting). In general, any feasible
 *   point suffices, but it might happen that the separation fails if the point is not in the relative interior.
 * - If @p primal is supplied the procedure tries to find a cut that is satisfied with equality by this point.
 *
 * @pre The trivially lifted cuts are valid! Otherwise, one can specify that the cuts are only locally valid (they would have to be lifted otherwise).
 *
 *
 * @b Description
 *
 * Consider some optimization problem over a polyhedron \f$P\f$ and some point \f$x^* \notin P\f$. If there is an
 * optimization oracle for this problem, valid cuts can be generated as follows.
 *
 * By the Minkowski-Weyl theorem we can write
 * \f[
 *     P = \mbox{conv}\{v \;:\; v \in V\} + \mbox{cone}\{r \;:\; r \in R\},
 * \f]
 * where \f$V\f$ is the (finite) set of vertices of \f$P\f$ and \f$R\f$ is the (finite) set of (extreme) rays of
 * \f$P\f$. Now consider any point \f$x^0 \in P\f$ and a primal separation point \f$p\f$. We can the consider the
 * following linear program:
 * \f[
 *     \begin{array}{lrl}
 *       \min & \gamma\phantom{\displaystyle\sum_{v \in V} v\, \lambda_v + \displaystyle\sum_{r \in R} r\, \mu_r + (x^* - x^0)\, \gamma + p} \\[1ex]
 *            & \displaystyle\sum_{v \in V} v\, \lambda_v + \displaystyle\sum_{r \in R} r\, \mu_r + (x^* - x^0)\, \gamma + p\, \delta & = x^*,\\[1ex]
 *            & - \displaystyle\sum_{v \in V} \lambda_v - \delta & = -1,\\[2ex]
 *            & \lambda, \mu, \gamma \geq 0,\; \delta \mbox{ free}. & \\
 *     \end{array}
 * \f]
 * Using \f$a_0\f$ as the dual variable for the last equation and \f$a\f$ for the other constraints, the dual problem is:
 * \f[
 *     \begin{array}{lr@{\;}ll}
 *       \max & (x^*)^T a - a_0 &&\\[1ex]
 *            & v^T a & \leq a_0 & \forall\, v \in V,\\[1ex]
 *            & r^T a & \leq 0   & \forall\, r \in R,\\[1ex]
 *            & (x^* - x^0)^T a & \leq 1, & \\[1ex]
 *            & p^T a & = a_0, & \\[1ex]
 *            & a, a_0 \mbox{ free}. & \\
 *     \end{array}
 * \f]
 * We can now interpret these two LPs and start with the dual:
 * - The dual variables \f$a\f$ and \f$a_0\f$ define an inequality \f$a^T x \leq a_0\f$ that is valid for \f$P\f$, because
 *   it is satisfied by all vertices and rays of \f$P\f$.
 * - The dual LP maximizes the violation \f$a_0 - a^T x^*\f$ of the inequality.
 * - The inequality \f$a^T x \leq a_0\f$ is satisfied with equality by the primal separation point \f$p\f$, thus
 *   providing primal separation.
 * - The normalizing inequality \f$(x^* - x^0)^T a \leq 1\f$ is used to bound the coefficients of the inequality.
 *
 * For the primal we observe the following:
 * - If for the optimal solution \f$\gamma = 0\f$ (and \f$\delta = 0\f$), then \f$x^*\f$ can be written as the convex and conic combination of vertices
 *   and rays of \f$P\f$. Thus, \f$x^* \in P\f$ and the point cannot be separated.
 * - If for the optimal solution \f$\gamma = 1\f$ (and \f$\delta = 0\f$), then \f$x^0\f$ can be written as the convex and conic combination of vertices
 *   and rays of \f$P\f$. Thus, we have proved that \f$x^0 \in P\f$. In this case \f$x^0\f$ should be chosen differently.
 *   (The line containing \f$x^0\f$ and \f$x^*\f$ does not contain points in the relative interior of
 *   \f$P\f$. Thus, if \f$x^0\f$ is in the relative interior and \f$x^*\f$ lies in the affine hull of \f$P\f$,
 *   this situation cannot occur.)
 *
 * Consequently, if the optimal \f$\gamma\f$ satisfies \f$0 < \gamma < 1\f$, this provides a violated
 * inequality. Moreover, this inequality provides a facet, if \f$(a, a_0)\f$ is a basic solution and we do not use a
 * primal separation point. This is holds since in this case the number of inequalities in the primal that are satisfied
 * with equality is equal to the dimension. Thus, the corresponding vertices and rays are affinely independent.
 *
 * <b> Solving the Primal </b>
 *
 * The primal problem can be solved in two ways: Either all vertices and rays can be enumerated or they can be produced
 * dynamically, i.e., in a column generation procedure.
 *
 * The pricing problem in the column generation then is: Find a vertex \f$v\f$ or ray \f$r\f$ of \f$P\f$ such that the
 * corresponding inequalities \f$v^T a \leq a_0\f$ and \f$r^T a \leq 0\f$ are violated. This amounts to maximize \f$a\f$
 * over \f$P\f$. If the result is larger than \f$a_0\f$ (or 0 in the case of rays), we have found a new column to be
 * added to the master problem. In the other case, we take the dual solution and produce a violated cut (or decide that
 * \f$x^* \in P\f$ or we have to change \f$x^0\f$).
 *
 * In terms of reduced costs the previous paragraph reads as follows: Let \f$y = (a, a_0)\f$ be the dual variables
 * w.r.t. the primal problem. If we denote the objective by \f$c\f$, the variables by \f$x\f$, and the system by \f$Ax =
 * b,\; x \geq 0\f$, Then the reduced cost of a variable \f$\lambda_v\f$ (corresponding to vertex \f$v \in V\f$):
 * \f[
 *     c_j - y^T A_j = 0 - a^T\, v + a_0.
 * \f]
 * Since we are minimizing, the reduced cost should be negative for a variable to be interesting. Thus \f$- a^T\, v +
 * a_0 < 0\f$ or equivalently \f$a^T\, v - a_0 > 0\f$. If the dual is unbounded (i.e., the primal unbounded), we get a
 * dual ray \f$(a, a_0)\f$. Since the objective coefficient of the interesting variables if 0, the pricing problem in
 * this case amounts to exactly the same problem as for the feasible case, except that we use a dual ray instead of a
 * dual solution.
 *
 * The key fact needed is that we need to solve a linear objective over the original problem, i.e., in general, the
 * problem we wanted to solve from the beginning. The idea in this context is to take a small subproblem of the original
 * and apply the above procedure to the subproblem. If the subproblem is small enough it can be solved in a practically
 * acceptable time. Thus, we separate cuts by optimizing, as the Ellipsoid Algorithm.
 *
 * Note that the produced cuts either need to be trivially liftable or be lifted by secondary procedure. In the
 * implementation <b>we assume trivial liftability</b>!
 *
 * <b> Local Cuts </b>
 *
 * The above idea first appeared as "local cuts" in the literature (Applegate, Bixby, Chvatal, Cook 2006). Their
 * formulation is somewhat different from the above:
 * \f[
 *     \begin{array}{lrl}
 *       \min & \displaystyle\sum_{i=1}^d (\sigma^+_i + \sigma^-_i)
 *              \phantom{\displaystyle\sum_{v \in V} v\, \lambda_v + \displaystyle\sum_{r \in R} r\, \mu_r + p\, \delta} \\[1ex]
 *            & \displaystyle\sum_{v \in V} v\, \lambda_v + \displaystyle\sum_{r \in R} r\, \mu_r + p\, \delta + \sigma^+ - \sigma^- & = x^*\\[1ex]
 *            & - \displaystyle\sum_{v \in V} \lambda_v - \delta & = -1\\[2ex]
 *            & \lambda, \mu, \sigma^+, \sigma^- \geq 0,\; \delta \mbox{ free} & \\
 *     \end{array}
 * \f]
 * If the objective value is 0, \f$x^*\f$ lies in the convex hull of the vertices and rays of \f$P\f$. Otherwise, we get
 * a valid cut, which not necessarily defines a facet. If the parameter @p primal is supplied and \f$p\f$ denotes the
 * separation point, then when the objective function is 0 the point \f$x^*\f$ lies in the convex hull of the vertices
 * and rays of \f$P\f$ and if \f$\delta>0\f$ the vertices of such a convex hull always inlcude \f$p\f$.
 *
 * Using as above \f$a_0\f$ as the dual variable for the last equation and \f$a\f$ for the other constraints, the dual problem is:
 * \f[
 *     \begin{array}{lr@{\;}ll}
 *       \max & (x^*)^T a - a_0 &&\\[1ex]
 *            & v^T a & \leq a_0 & \forall\, v \in V\\[1ex]
 *            & r^T a & \leq 0   & \forall\, r \in R\\[1ex]
 *            & p^T a & = a_0 & \\[1ex]
 *            & -1 \leq a_i & \leq 1 & \forall\, i = 1, \dots, d\\[1ex]
 *            & a, a_0 \mbox{ free} & \\
 *     \end{array}
 * \f]
 *
 * Analogously, if the optimal value of the dual problem is positive, then the inequality \f$ax\leq a_0\f$ is valid
 * (since it is satisfied by all vertices of the polyhdron, has nonpositive inner product with all its rays and, in case
 * the parameter @p primal is supplied, is satisfied by the point \f$p\f$ at equality. Moreover such an inequality is
 * violated by the point \f$x^*\f$ by an amount that is given by the optimal objective function value.
 *
 * <b>Implementation Details</b>
 *
 * In the implementation takes care of the following issues:
 * - As explained above, the primal master problem might be infeasible in between. In this case we perform Farkas
 *   pricing, i.e., use a dual ray for the pricing problem.
 * - We try to make "nice" cuts out of the resulting dual variables. This is done via a rational approximation. The cuts
 *   are directly added to SCIP.
 */
SCIP_RETCODE LC_generateCut(
   SCIP*                 scip,               /**< SCIP pointer */
   int                   dim,                /**< dimension of subproblem */
   const SCIP_Real*      x,                  /**< point to be separated in reduced space */
   SCIP_Real*            xrelint,            /**< relative interior point (or NULL) */
   SCIP_VAR**            vars,               /**< original variables corresponding to subproblem */
   LC_ORACLE((*Oracle)),                     /**< optimization oracle callback */
   const LC_DATA*        data,               /**< data passed to oracle (or NULL) */
   int                   nsols,              /**< number of initial feasible solutions */
   const SCIP_Real**     sols,               /**< initial feasible solutions */
   const SCIP_Real*      primal,             /**< primal separation point or NULL if primal separation is turned off */
   int                   ninitialrows,       /**< number of initial rows to be added to master LP */
   SCIP_Real*            initialrowslhs,     /**< left hand sides of initial rows */
   SCIP_Real*            initialrowsrhs,     /**< right hand sides of initial rows */
   int*                  initialrowsbeg,     /**< beginning indices of initial rows */
   SCIP_Real*            initialrowsvals,    /**< values of initial rows */
   int*                  initialrowsinds,    /**< variable indices of initial rows */
   SCIP_Bool             targetcut,          /**< whether the separation is based on a target cuts (local cuts otherwise) */
   SCIP_Bool             tilting,            /**< whether inequality tilting is performed (only for local cuts) */
   SCIP_Bool             checksols,          /**< whether the solutions should be checked for feasibility in the orginal problem */
   SCIP_Bool             nicelocalcuts,      /**< try to make a nice cut by finding a rational approximation */
   SCIP_Real*            cutcoefs,           /**< cut coefficients of (tilted) cut (or NULL) */
   SCIP_Real*            cutrhs,             /**< pointer to store the rhs of the cut (or NULL) */
   SCIP_Bool*            success,            /**< pointer to store whether cut generation was successful */
   const SCIP_Real*      maxtime,            /**< maximal time that should be used for separation (or NULL) */
   SCIP_Real*            usedtime            /**< returns the used time */
   )
{
   /* basic data */
   SCIP_LPI* masterlp;                       /* the lp used to solve the master problems (via LPI) */
   SCIP_Real starttime;                      /* start time */
   SCIP_CLOCK* timeLP;                       /* clock for measuring the time for solving LPs */
   SCIP_CLOCK* timeOracle;                   /* clock for measuring the time for calling the oracle */
#ifdef SCIP_DEBUG
   int nmasterlps = 0;                       /* number of master LPs solved */
#endif
   int nmastercols = 0;                      /* number of columns in master LP */
   SCIPdebug( int noldmastercols; )

   SCIP_Real colobj;                         /* value of the colum objective used for pricing */
   SCIP_Real* relint;                        /* relative interior point */
   SCIP_Real timelimit;
   SCIP_Real val;
   int lpdim;
   int cnt;
   int i;
   int k;
   int j;

   /* data for communication with oracle */
   SCIP_Real* oracleobj;                     /* objective function for oracle */
   SCIP_Real* dualsol;                       /* dual solution of master LP */
   int noraclesols;                          /* number of oracle solutions returned */
   int nmaxoraclesols = MAXNSOLS;            /* maximal number of oracle solutions that can be returned */
   SCIP_Real** oraclesols;                   /* storage for solutions of oracle */
   SCIP_Real* oracleobjvals;                 /* objective function values of solutions of oracle */

   /* storage for created points */
   SCIP_Bool storepoints = FALSE;            /* whether primal points should be stored */
   int npoints = 0;                          /* number of created points */
   int nmaxpoints = 0;                       /* maximal number of points that can be stored */
   int pointsize = 0;                        /* number of nonzeros of points stored */
   int maxpointsize = 0;                     /* maximal number of nonzeros that can be stored */
   int* pointbeg = NULL;                     /* index of points in pointind/pointval arrays */
   int* pointind = NULL;                     /* indices of nonzero components of created points */
   SCIP_Real* pointval = NULL;               /* values of nonzero components of created points */

   /* storage for creating columns */
   SCIP_Real obj;
   SCIP_Real lb;
   SCIP_Real ub;
   int matbeg;
   int* matind;
   SCIP_Real* matval;
   SCIP_Real* rhs = NULL;
   SCIP_Real* lhs = NULL;
   SCIP_Real masterObjval;
   SCIP_Bool perturbed = FALSE;

   assert( scip != NULL );
   assert( dim > 0 );
   assert( x != NULL );
   assert( vars != NULL );
   assert( Oracle != NULL );
   assert( usedtime != NULL );
   assert( success != NULL );
   assert( nsols == 0 || sols != NULL );

   assert( ninitialrows == 0 || initialrowslhs != NULL );
   assert( ninitialrows == 0 || initialrowsrhs != NULL );
   assert( ninitialrows == 0 || initialrowsbeg != NULL );
   assert( ninitialrows == 0 || initialrowsvals != NULL );
   assert( ninitialrows == 0 || initialrowsinds != NULL );

   SCIPdebugMsg(scip, "Handle local cut separation with dim=%d ...\n", dim);

   /* initializations */
   *usedtime = 0.0;
   *success = FALSE;

   /* check whether relative interior point has been supplied for target cuts */
   if ( targetcut && xrelint == NULL && sols == NULL )
   {
      SCIPerrorMessage("Cannot call target cuts without relative interior point or primal solutions to compute it.\n");
      return SCIP_INVALIDDATA;
   }

   /* init timelimit */
   starttime = SCIPgetSolvingTime(scip);
   SCIP_CALL( SCIPgetRealParam(scip, "limits/time", &timelimit) );
   if ( ! SCIPisInfinity(scip, timelimit) )
      timelimit -= SCIPgetSolvingTime(scip);
   if ( maxtime != NULL && *maxtime < timelimit )
      timelimit = MAX(*maxtime, 0.0);

   /* create clocks */
   SCIP_CALL( SCIPcreateClock(scip, &timeLP) );
   SCIP_CALL( SCIPcreateClock(scip, &timeOracle) );

   /* get space for points */
   if ( tilting && xrelint == NULL )
   {
      storepoints = TRUE;

      nmaxpoints = INIT_NMAXPOINTS;
      assert( nmaxpoints > 0 );
      maxpointsize = nmaxpoints * dim;

      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &pointbeg, nmaxpoints) );
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &pointind, maxpointsize) );
      SCIP_CALL( SCIPallocBlockMemoryArray(scip, &pointval, maxpointsize) );
   }

   /* allocate space */
   if ( targetcut )
      lpdim = dim;
   else
      lpdim = dim + 1;

   SCIP_CALL( SCIPallocBufferArray(scip, &dualsol, lpdim + ninitialrows) );
   SCIP_CALL( SCIPallocBufferArray(scip, &oracleobj, dim) );
   SCIP_CALL( SCIPallocBufferArray(scip, &matind, lpdim + ninitialrows) );
   SCIP_CALL( SCIPallocBufferArray(scip, &matval, lpdim + ninitialrows) );

   /* possibly initialize relative interior point */
   if ( xrelint == NULL )
   {
      SCIP_CALL( SCIPallocBufferArray(scip, &relint, dim) );

      for (j = 0; j < dim; ++j)
         relint[j] = 0.0;

      /* if there are enough initial points, compute their barycenter */
      if ( nsols > 0 )
      {
         for (i = 0; i < nsols; ++i)
         {
            for (j = 0; j < dim; ++j)
               relint[j] += sols[i][j];
         }
         for (j = 0; j < dim; ++j)
            relint[j] = relint[j] / ((SCIP_Real) nsols);
      }
   }
   else
      relint = xrelint;
   assert( relint != NULL );

   /* create master problem */
   SCIP_CALL( SCIPlpiCreate(&masterlp, SCIPgetMessagehdlr(scip), "master", SCIP_OBJSEN_MINIMIZE) );

   /* prepare rows */
   SCIP_CALL( SCIPallocBufferArray(scip, &lhs, lpdim) );
   SCIP_CALL( SCIPallocBufferArray(scip, &rhs, lpdim) );
   for (i = 0; i < dim; ++i)
   {
      if ( targetcut )
      {
         assert( relint != NULL );
         lhs[i] = x[i] - relint[i];
         rhs[i] = lhs[i];
      }
      else
      {
         lhs[i] = x[i];
         rhs[i] = x[i];
      }
   }

   if ( ! targetcut )
   {
      assert( lpdim == dim + 1 );
      lhs[dim] = -1.0;
      rhs[dim] = -1.0;
   }

   /* create rows */
   SCIP_CALL( SCIPlpiAddRows(masterlp, lpdim, lhs, rhs, NULL, 0, NULL, NULL, NULL) );

   /* add rows for initial rows */
   if ( ninitialrows > 0 )
   {
      SCIP_CALL( SCIPlpiAddRows(masterlp, ninitialrows, initialrowslhs, initialrowsrhs, NULL, 0, NULL, NULL, NULL) );
   }

   /* init parameters for columns */
   matbeg = 0;

   /* add primal separation column if required */
   if ( primal != NULL )
   {
      cnt = 0;

      obj = targetcut ? 1.0 : 0.0;
      lb = -SCIPlpiInfinity(masterlp);
      ub = SCIPlpiInfinity(masterlp);

      for (j = 0; j < dim; ++j)
      {
         val = primal[j];
	 if ( ! SCIPisFeasZero(scip, val) )
	 {
	    matind[cnt] = j;
	    matval[cnt++] = val;
	 }
      }
      if ( ! targetcut )
      {
         matind[cnt] = dim;
         matval[cnt++] = -1.0;
      }

      SCIP_CALL( SCIPlpiAddCols(masterlp, 1, &obj, &lb, &ub, NULL, cnt, &matbeg, matind, matval) );
      ++nmastercols;
   }

   /* add initial solutions */
   if ( nsols > 0 )
   {
      obj = targetcut ? 1.0 : 0.0;
      lb = 0.0;
      ub = SCIPlpiInfinity(masterlp);

      for (k = 0; k < nsols; ++k)
      {
         cnt = 0;
         for (j = 0; j < dim; ++j)
         {
            if ( targetcut )
               val = sols[k][j] - relint[j];
            else
               val = sols[k][j];

            if ( ! SCIPisFeasZero(scip, val) )
            {
               matind[cnt] = j;
               matval[cnt++] = val;
            }
         }

         if ( ! targetcut )
         {
            matind[cnt] = dim;
            matval[cnt++] = -1.0;
         }

         /* possibly add initial row part */
         for (i = 0; i < ninitialrows; ++i)
         {
#ifndef NDEBUG
            SCIP_Bool isintegral = TRUE;
#endif

            val = 0.0;
            for (j = initialrowsbeg[i]; j < initialrowsbeg[i+1]; ++j)
            {
               int ind;
               ind = initialrowsinds[j];
               assert( 0 <= ind && ind < dim );

               val += initialrowsvals[j] * sols[k][ind];

#ifndef NDEBUG
               /* check integrality */
               if ( SCIPvarIsIntegral(vars[ind]) )
               {
                  assert( SCIPisFeasIntegral(scip, sols[k][ind]) );
                  if ( ! SCIPisIntegral(scip, initialrowsvals[j]) )
                     isintegral = FALSE;
               }
               else
                  isintegral = FALSE;
#endif
            }
            assert( ! isintegral || SCIPisFeasIntegral(scip, val) );

            if ( ! SCIPisFeasZero(scip, val) )
            {
               matind[cnt] = lpdim + i;
               matval[cnt++] = val;
            }
         }

         SCIP_CALL( SCIPlpiAddCols(masterlp, 1, &obj, &lb, &ub, NULL, cnt, &matbeg, matind, matval) );
         ++nmastercols;
      }
      SCIPdebugMsg(scip, "Added initial columns: %d\n", nsols);
   }

   /* add 2 * dim variables for norm L1 minimization */
   if ( ! targetcut )
   {
      obj = 1.0;
      lb = 0.0;
      ub = SCIPlpiInfinity(masterlp);

      matval[0] = 1.0;
      for (i = 0; i < dim; ++i)
      {
         matind[0] = i;
         SCIP_CALL( SCIPlpiAddCols(masterlp, 1, &obj, &lb, &ub, NULL, 1, &matbeg, matind, matval) );
         ++nmastercols;
      }
      matval[0] = -1.0;
      for (i = 0; i < dim; ++i)
      {
         matind[0] = i;
         SCIP_CALL( SCIPlpiAddCols(masterlp, 1, &obj, &lb, &ub, NULL, 1, &matbeg, matind, matval) );
         ++nmastercols;
      }
   }

#ifdef SCIP_CHECK
   SCIP_CALL( SCIPlpiWriteLP(masterlp, "starting_master.lp") );
#endif

   /* for debugging */
#ifdef SCIP_DEBUG
   for (j = 0; j < lpdim + ninitialrows; ++j)
      dualsol[j] = -1.0;
   for (j = 0; j < dim; ++j)
      oracleobj[j] = -1.0;
#endif

   /* prepare space for additional solutions */
   SCIP_CALL( SCIPallocBufferArray(scip, &oraclesols, nmaxoraclesols) );
   SCIP_CALL( SCIPallocBufferArray(scip, &oracleobjvals, nmaxoraclesols) );
   for (i = 0; i < nmaxoraclesols; ++i)
   {
      SCIP_CALL( SCIPallocBufferArray(scip, &(oraclesols[i]), dim) );
   }

   /* prepare settings for new columns */
   obj = targetcut ? 1.0 : 0.0;
   lb = 0.0;
   ub = SCIPlpiInfinity(masterlp);

   /* main tilting loop */
   do
   {
      SCIP_Bool pricedcol = FALSE;              /* Has a column be priced out? */

      if ( tilting )
      {
         SCIPdebugMsg(scip, "tilting loop ...\n");
         perturbed = FALSE;
      }

      if ( SCIPisStopped(scip) )
         goto terminate;

      /* main loop */
      do
      {
         SCIP_Real t;
         SCIP_Bool farkaspricing = FALSE;

         t = SCIPgetSolvingTime(scip) - starttime;
         pricedcol = FALSE;

         /* check whether we reach the time limit */
         if ( t > timelimit )
            goto terminate;

         /* set time limit */
         if ( ! SCIPisInfinity(scip, timelimit) )
         {
            SCIP_CALL( SCIPlpiSetRealpar(masterlp, SCIP_LPPAR_LPTILIM, timelimit - t) );
            /* SCIPdebugMsg(scip, "Solving LP %u with time limit %g\n", nmasterlps, timelimit - t); */
         }
         /* else
            SCIPdebugMsg(scip, "Solving LP %u.\n", nmasterlps); */

#ifndef NDEBUG
         {
            int nRows, nCols;
            SCIP_CALL( SCIPlpiGetNRows(masterlp, &nRows) );
            SCIP_CALL( SCIPlpiGetNCols(masterlp, &nCols) );
            assert( nRows == lpdim + ninitialrows );
            assert( nCols == nmastercols );
         }
#endif
#ifdef SCIP_CHECK
         SCIP_CALL( SCIPlpiWriteLP(masterlp, "master.lp") );
#endif

         /* optimize master problem (dual is usually faster than primal, even for column generation) */
         SCIP_CALL( SCIPstartClock(scip, timeLP) );
         SCIP_CALL( SCIPlpiSolveDual(masterlp) );
         SCIP_CALL( SCIPstopClock(scip, timeLP) );

         /* check whether we reached the time limit */
         if ( SCIPlpiIsTimelimExc(masterlp) )
            goto terminate;

         /* in primal separation the LP might become unbounded (the dual is infeasible) - in this case we terminate separation */
         if ( primal != NULL && SCIPlpiIsDualInfeasible(masterlp) )
         {
            SCIPdebugMsg(scip, "The dual is infeasible, using primal separation -> stop separation.\n");
            goto terminate;
         }

         /* check status and determine whether we have to perform Farkas pricing */
         if ( SCIPlpiIsOptimal(masterlp) )
         {
            /* get dual variables */
            SCIP_CALL( SCIPlpiGetSol(masterlp, &masterObjval, NULL, dualsol, NULL, NULL) );
         }
         else if ( SCIPlpiIsDualUnbounded(masterlp) || SCIPlpiIsPrimalInfeasible(masterlp) )
         {
            if ( ! SCIPlpiHasDualRay(masterlp) )
            {
               SCIPerrorMessage("solver cannot return dual ray - have to exit (status: %d)\n", SCIPlpiGetInternalStatus(masterlp) );
               return SCIP_ERROR;
            }

            /* get dual ray */
            SCIP_CALL( SCIPlpiGetDualfarkas(masterlp, dualsol) );
            masterObjval = SCIPlpiInfinity(masterlp);
            farkaspricing = TRUE;
         }
         else
         {
            SCIP_CALL( SCIPlpiWriteLP(masterlp, "err.lp") );
            SCIPerrorMessage("error in solving LP (status: %d)\n", SCIPlpiGetInternalStatus(masterlp) );
            return SCIP_ERROR;
         }
#ifdef SCIP_DEBUG
         ++nmasterlps;
#endif

         /* note: at this time dualsol contains the dual values or a dual ray */

         /* check for violated cuts: positive objective means that we cannot represent the LP-relaxation solution, i.e., can separate a cut */
         colobj = targetcut ? 1.0 : 0.0;
         if ( SCIPisEfficacious(scip, masterObjval - colobj) )
         {
            LC_STATUS status;                       /* Status of oracle */

            /* determine rhs */
            if ( targetcut )
            {
               if ( farkaspricing )
                  colobj = 0.0;
               else
                  assert( colobj == 1.0 );

               for (j = 0; j < dim; ++j)
                  colobj += dualsol[j] * relint[j];
            }
            else
               colobj = dualsol[dim];

            /* set up objective for oracle */
            for (j = 0; j < dim; ++j)
               oracleobj[j] = dualsol[j];

            /* add initial row contribution */
            for (i = 0; i < ninitialrows; ++i)
            {
#if 0
               /* The dual variables can have any sign for equations. They should be nonnegative for >= inequalities and nonpositive for <= inequalities */
               assert( SCIPisFeasGE(scip, dualsol[lpdim + i], 0.0) || SCIPisEQ(scip, initialrowslhs[i], initialrowsrhs[i]) || SCIPisInfinity(scip, -initialrowslhs[i]) );
               assert( SCIPisFeasLE(scip, dualsol[lpdim + i], 0.0) || SCIPisEQ(scip, initialrowslhs[i], initialrowsrhs[i]) || SCIPisInfinity(scip, initialrowsrhs[i]) );
#endif
               for (j = initialrowsbeg[i]; j < initialrowsbeg[i+1]; ++j)
               {
                  assert( 0 <= initialrowsinds[j] && initialrowsinds[j] < dim );
                  oracleobj[initialrowsinds[j]] += initialrowsvals[j] * dualsol[lpdim + i];
               }
            }

            /* call oracle */
            SCIP_CALL( SCIPstartClock(scip, timeOracle) );
            SCIP_CALL( Oracle(scip, dim, data, oracleobj, colobj, nmaxoraclesols, &noraclesols, oraclesols, oracleobjvals, &status) );
            SCIP_CALL( SCIPstopClock(scip, timeOracle) );

            if ( status != LC_STATUS_OPTIMAL && status != LC_STATUS_CUTOFF )
            {
               SCIPerrorMessage("Something went wrong in the oracle (status: %d).\n", status);
               goto terminate;
            }
            /* note: noraclesols == 0 is possible if the cutoff is set */
            assert( noraclesols <= nmaxoraclesols );
            assert( noraclesols > 0 || status == LC_STATUS_CUTOFF );

            /* add solutions if they are worth it */
            SCIPdebug( noldmastercols = nmastercols; )
            for (k = 0; k < noraclesols; ++k)
            {
               if ( SCIPisEfficacious(scip, oracleobjvals[k] - colobj) )
               {
#ifndef NDEBUG
                  /* check whether computed solution value was correct */
                  {
                     val = 0.0;
                     for (j = 0; j < dim; ++j)
                        val += oracleobj[j] * oraclesols[k][j];
                     assert( SCIPisFeasEQ(scip, val, oracleobjvals[k]) );
                  }
#endif

                  /* try whether solution is feasible for original problem */
                  if ( checksols )
                  {
                     SCIP_CALL( checkPoint(scip, oraclesols[k], dim, vars) );
                  }

                  /* prepare new column */
                  cnt = 0;
                  for (j = 0; j < dim; ++j)
                  {
                     if ( targetcut )
                        val = oraclesols[k][j] - relint[j];
                     else
                        val = oraclesols[k][j];

                     if ( ! SCIPisZero(scip, val) )
                     {
                        matind[cnt] = j;
                        matval[cnt++] = val;
                     }
                  }

                  if ( ! targetcut )
                  {
                     matind[cnt] = dim;
                     matval[cnt++] = -1.0;
                  }
                  pricedcol = TRUE;

                  /* possibly add initial row part */
                  for (i = 0; i < ninitialrows; ++i)
                  {
#ifndef NDEBUG
                     SCIP_Bool isintegral = TRUE;
#endif

                     val = 0.0;
                     for (j = initialrowsbeg[i]; j < initialrowsbeg[i+1]; ++j)
                     {
                        int ind;
                        ind = initialrowsinds[j];
                        assert( 0 <= ind && ind < dim );

                        val += initialrowsvals[j] * oraclesols[k][ind];

#ifndef NDEBUG
                        /* check integrality */
                        if ( SCIPvarIsIntegral(vars[ind]) )
                        {
                           assert( SCIPisFeasIntegral(scip, oraclesols[k][ind]) );
                           if ( ! SCIPisIntegral(scip, initialrowsvals[j]) )
                              isintegral = FALSE;
                        }
                        else
                           isintegral = FALSE;
#endif
                     }
                     assert( ! isintegral || SCIPisFeasIntegral(scip, val) );

                     if ( ! SCIPisZero(scip, val) )
                     {
                        matind[cnt] = lpdim + i;
                        matval[cnt++] = val;
                     }
                  }

                  /* add column */
                  SCIP_CALL( SCIPlpiAddCols(masterlp, 1, &obj, &lb, &ub, NULL, cnt, &matbeg, matind, matval) );
                  ++nmastercols;
               }
            }

            /* possibly store points */
            if ( storepoints && noraclesols > 0 )
            {
               assert( pointbeg != NULL );
               assert( pointind != NULL );
               assert( pointval != NULL );

               /* ensure space */
               SCIP_CALL( ensurePointNumber(scip, &pointbeg, &nmaxpoints, npoints + noraclesols + 1) );
               SCIP_CALL( ensurePointSize(scip, &pointind, &pointval, &maxpointsize, pointsize + noraclesols * dim) );

               for (k = 0; k < noraclesols; ++k)
               {
                  if ( SCIPisEfficacious(scip, oracleobjvals[k] - colobj) )
                  {
                     pointbeg[npoints++] = pointsize;

#ifdef SCIP_CHECK
                     checkPointIsNew(scip, oraclesols[k], dim, npoints-1, pointbeg, pointind, pointval);
#endif

                     /* store solution as new point */
                     for (j = 0; j < dim; ++j)
                     {
                        val = oraclesols[k][j];
                        if ( ! SCIPisFeasZero(scip, val) )
                        {
                           pointind[pointsize] = j;
                           pointval[pointsize++] = val;
                        }
                     }
                     assert( pointsize < maxpointsize );
                  }
               }
               /* set pointbeg one beyond the last (space reserved above) */
               pointbeg[npoints] = pointsize;

               assert( npoints <= nmastercols );
               assert( npoints < nmaxpoints );
               assert( pointsize < maxpointsize );
            }

#ifdef SCIP_DEBUG
            if ( noraclesols == 0 )
               SCIPinfoMessage(scip, NULL, "The oracle found no column to be priced.\n");
            else
            {
               if ( ! pricedcol )
                  SCIPinfoMessage(scip, NULL, "No solution of the oracle creates a column to be priced.\n");
            }

            if ( SCIPisEfficacious(scip, oracleobjvals[0] - colobj) && (nmasterlps % OUTPUTFREQ == 0) )
            {
               SCIPinfoMessage(scip, NULL, "iter: %4d, slack: %8.6g, #cols: %4d (colobj: %8.5g, added columns: %2d; best pricing violation = %f)\n",
                  nmasterlps, masterObjval, nmastercols, colobj, nmastercols - noldmastercols, REALABS(colobj-oracleobjvals[0]));
            }
#endif
         }
      }
      while ( pricedcol );

#ifdef SCIP_DEBUG
      /* final output */
      SCIPinfoMessage(scip, NULL, "final iter: %4d, obj: %8.6g, #cols: %4d\n", nmasterlps, masterObjval, nmastercols);

      if ( primal != NULL )
      {
         SCIP_Real activity = 0.0;

         for (j = 0; j < dim; ++j)
         {
            if ( targetcut )
               activity += (primal[j] - relint[j]) * dualsol[j];
            else
               activity += primal[j] * dualsol[j];
         }
         assert( SCIPisFeasEQ(scip, activity, targetcut ? 1.0 : dualsol[dim]) );
      }
#endif

      /* check if the solution obtained provides a valid local cut (infeasible master LPs always lead to cut) */
      colobj = targetcut ? 1.0 : 0.0;
      if ( SCIPisEfficacious(scip, masterObjval - colobj) )
      {
         SCIP_CALL( LC_checkCut(scip, nicelocalcuts, x, dim, dualsol, MAXINTEGER, MAXDENOMINATOR, MAXABSCOEF, Oracle, data,
               oraclesols, oracleobjvals, cutcoefs, cutrhs, success) );

         /* possibly call tilting */
         if ( tilting && ! targetcut )
         {
            if ( xrelint == NULL )
            {
               assert( storepoints );
               assert( pointbeg != NULL );
               assert( pointind != NULL );
               assert( pointval != NULL );

               /* recompute relative interior point and change lp */
               SCIPdebugMsg(scip, "Recomputing relative interior point ...\n");

               /* initialize relative interior point */
               for (j = 0; j < dim; ++j)
                  relint[j] = 0.0;

               /* compute barycenter */
               assert( 0 < npoints && npoints < nmaxpoints );
               assert( pointbeg[npoints] == pointsize );
               for (i = 0; i < npoints; ++i)
               {
                  for (j = pointbeg[i]; j < pointbeg[i+1]; ++j)
                  {
                     assert( 0 <= pointind[j] && pointind[j] < dim );
                     relint[pointind[j]] += pointval[j];
                  }
               }
               /* divide */
               assert( npoints > 0 );
               for (j = 0; j < dim; ++j)
                  relint[j] = relint[j] / ((SCIP_Real) npoints);
            }

            SCIP_CALL( perturb(scip, masterlp, dim, x, relint, dualsol, &perturbed) );
            if ( perturbed )
               SCIPdebugMsg(scip, "Perturbed problem for tilting ...\n");
         }
      }
   }
   while ( perturbed );

 terminate:
   for (k = 0; k < nmaxoraclesols; ++k)
      SCIPfreeBufferArrayNull(scip, &oraclesols[k]);
   SCIPfreeBufferArrayNull(scip, &oracleobjvals);
   SCIPfreeBufferArrayNull(scip, &oraclesols);

   SCIPfreeBufferArray(scip, &rhs);
   SCIPfreeBufferArray(scip, &lhs);

   if ( xrelint == NULL )
      SCIPfreeBufferArray(scip, &relint);

   SCIPfreeBufferArray(scip, &matval);
   SCIPfreeBufferArray(scip, &matind);
   SCIPfreeBufferArray(scip, &oracleobj);
   SCIPfreeBufferArray(scip, &dualsol);

   SCIP_CALL( SCIPlpiFree(&masterlp) );

   if ( storepoints )
   {
      SCIPfreeBlockMemoryArray(scip, &pointbeg, nmaxpoints);
      SCIPfreeBlockMemoryArray(scip, &pointind, maxpointsize);
      SCIPfreeBlockMemoryArray(scip, &pointval, maxpointsize);
   }
   else
   {
      assert( nmaxpoints == 0 );
      assert( maxpointsize == 0 );
      assert( pointbeg == NULL );
      assert( pointind == NULL );
      assert( pointval == NULL );
   }

   *usedtime = SCIPgetTotalTime(scip) - starttime;

   SCIP_CALL( SCIPfreeClock(scip, &timeLP) );
   SCIP_CALL( SCIPfreeClock(scip, &timeOracle) );

   return SCIP_OKAY;
}
