infinity: Implement cone_bound_find_lambda

CMakeLists.txt

@@ -1,16 +1,18 @@
 cmake_minimum_required(VERSION 2.8)
 project(multirow)
 
-set(CMAKE_CXX_STANDARD 11)
-set(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Werror")
-set(CMAKE_C_FLAGS  "${CMAKE_CXX_FLAGS} -Werror")
-
 list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/cmake")
 find_package(CPLEX REQUIRED)
 include_directories(${CPLEX_INCLUDE_DIR})
 
 find_package(GMP REQUIRED)
 find_package(OpenMP REQUIRED)
+find_package(BLAS REQUIRED)
+find_package(LAPACKE REQUIRED)
+
+set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Werror ${OpenMP_CXX_FLAGS} -O3")
+set(CMAKE_C_FLAGS  "${CMAKE_CXX_FLAGS} -Werror ${OpenMP_C_FLAGS} -O3")
 
 include_directories(${gtest_SOURCE_DIR}/include)
 include_directories(infinity/library/include)

cmake/FindCPLEX.cmake

@@ -1,4 +1,6 @@
 set(CMAKE_LIBRARY_PATH ${CMAKE_LIBRARY_PATH}
+  $ENV{HOME}/lib/
+  /opt/cplex-12.6/cplex/lib/x86-64_linux/static_pic
   /software/cplex-12.6/distribution/cplex/lib/x86-64_linux/static_pic
   /Users/axavier/Applications/IBM/ILOG/CPLEX_Studio1262/cplex/lib/x86-64_osx/static_pic)
 
@@ -6,7 +8,9 @@ find_library(CPLEX_LIBRARIES
   NAMES cplex cplex1220 cplex1240 cplex1260 cplex1261 cplex1262)
 
 find_path(CPLEX_INCLUDE_DIR NAMES ilcplex/cplex.h PATHS
-  /Users/axavier/Applications/IBM/ILOG/CPLEX_Studio1262/cplex/include
+  $ENV{HOME}/include/
+  /opt/cplex-12.6/cplex/include/
+  /Users/axavier/Applications/IBM/ILOG/CPLEX_Studio1262/cplex/include/
   /software/cplex-12.6/distribution/cplex/include/)
 
 include(FindPackageHandleStandardArgs)

cmake/FindLAPACKE.cmake

@@ -0,0 +1,11 @@
+find_library(LAPACKE_LIBRARIES
+        NAMES lapacke)
+
+find_path(LAPACKE_INCLUDE_DIR NAMES lapacke.h)
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(LAPACKE
+        DEFAULT_MSG
+        LAPACKE_LIBRARIES
+        LAPACKE_INCLUDE_DIR)
+mark_as_advanced(LAPACKE_LIBRARIES LAPACKE_INCLUDE_DIR)

infinity/library/src/infinity-nd.c

@@ -22,6 +22,7 @@
 #include <multirow/util.h>
 
 #include <infinity/infinity-nd.h>
+#include <multirow/linalg.h>
 
 static long lp_count = 0;
 static double lp_time = 0;
@@ -41,6 +42,41 @@ static double sfree_mip_time = 0;
 static long scale_ahull_lp_count = 0;
 static double scale_ahull_lp_time = 0;
 
+static int cone_bound_find_lambda(const struct RayList *rays,
+                                  const double *f,
+                                  const double *x,
+                                  double *lambda)
+{
+    int rval = 0;
+    double *A = 0;
+    double *b = 0;
+    int dim = rays->dim;
+    int nrays = rays->nrays;
+
+    A = (double*) malloc(nrays * dim * sizeof(double));
+    b = (double*) malloc(dim * sizeof(double));
+    abort_if(!A, "could not allocate A");
+    abort_if(!b, "could not allocate b");
+
+    for (int i = 0; i < dim; i++)
+    {
+        b[i] = x[i] - f[i];
+        for (int j = 0; j < nrays; j++)
+            A[i * nrays + j] = rays->values[j * dim + i];
+    }
+
+    rval = LINALG_solve(nrays, dim, A, b, lambda);
+    abort_if(rval, "LINALG_solve failed");
+
+    for (int i = 0; i < dim; i++)
+        abort_iff(lambda[0] < 0, "lambda[i] is negative (%lf)", i, lambda[i]);
+
+CLEANUP:
+    if(A) free(A);
+    if(b) free(b);
+    return rval;
+}
+
 static int create_sfree_mip(int nrows,
                             int nrays,
                             const double *f,
@@ -365,10 +401,6 @@ static int create_tight_rays_lp(int nrows,
     rval = LP_relax(lp);
     abort_if(rval, "LP_relax failed");
 
-    //rval = LP_write(lp, "tight-rays.lp");
-    //abort_if(rval, "LP_write failed");
-
-
 CLEANUP:
     if(rmatind) free(rmatind);
     if(rmatval) free(rmatval);
@@ -430,9 +462,8 @@ static int create_violated_cone_lp(int nrows,
     rval = LP_relax(lp);
     abort_if(rval, "LP_relax failed");
 
-    //rval = LP_write(lp, "violated-cone.lp");
+//    rval = LP_write(lp, "violated-cone.lp");
-    //abort_if(rval, "LP_write failed");
+//    abort_if(rval, "LP_write failed");
-    //UTIL_pause();
 
 CLEANUP:
     if(rmatind) free(rmatind);
@@ -602,7 +633,7 @@ static int find_interior_point_enum(const int nrows,
                 }
             }
 
-    if(!DOUBLE_geq(best_value, 1)) *found = 1;
+    if(best_value < 0.999) *found = 1;
 
 CLEANUP:
     if(beta2) free(beta2);
@@ -635,14 +666,14 @@ static int find_interior_point_cplex(const int nrows,
     rval = create_sfree_mip(nrows, nrays, f, rays, beta, epsilon, &lp);
     abort_if(rval, "greate_sfree_mip failed");
 
-    log_debug("  solving sfree mip...\n");
+    log_debug("    solving sfree mip...\n");
     rval = LP_optimize(&lp, &infeasible);
     if(rval == ERR_MIP_TIMEOUT) goto CLEANUP;
     abort_if(rval, "LP_optimize failed");
 
     if(infeasible)
     {
-        log_debug("  mip is infeasible. Stopping.\n");
+        log_debug("        mip is infeasible. Stopping.\n");
         *found = 0;
         goto CLEANUP;
     }
@@ -663,11 +694,11 @@ static int find_interior_point_cplex(const int nrows,
     rval = LP_get_obj_val(&lp, &objval);
     abort_if(rval, "LP_get_obj_val failed");
 
-    log_debug("    obj = %.8lf\n", objval);
+    log_debug("        obj: %.8lf\n", objval);
 
-    if(DOUBLE_geq(objval, 1.0))
+    if(objval > 0.999)
     {
-        log_debug("  set is lattice-free\n");
+        log_debug("    set is lattice-free\n");
         *found = 0;
         goto CLEANUP;
     }
@@ -857,7 +888,7 @@ static int find_violated_cone(int nrows,
                               double *sbar,
                               int *violated_found)
 {
-    log_verbose("      find_violated_cone\n");
+    log_debug("    finding violated cone:\n");
     int rval = 0;
 
     struct LP lp;
@@ -894,30 +925,30 @@ static int find_violated_cone(int nrows,
     rval = LP_get_obj_val(&lp, &obj);
     abort_if(rval, "LP_get_obj_val failed");
 
-    log_verbose("  o=%.8lf\n", obj);
+    log_debug("      obj: %.8lf\n", obj);
 
-    if(DOUBLE_geq(obj, 0.999))
+    if(obj > 0.999)
     {
         *violated_found = 0;
+        log_debug("    no violated cone found\n");
     }
     else
     {
         *violated_found = 1;
 
-        log_verbose("Violated cone found\n");
+        log_debug("    violated cone found:\n");
-        log_verbose("  f=%.8lf %.8lf\n", f[0], f[1]);
+        log_debug("      f: %12.8lf %12.8lf\n", f[0], f[1]);
-        log_verbose("  x=%.8lf %.8lf\n", x[0], x[1]);
+        log_debug("      x: %12.8lf %12.8lf\n", x[0], x[1]);
 
         for(int i = 0; i < nrays; i++)
         {
             rx[i] = (sbar[i] > 1e-9);
 
-            if(rx[i]) if_verbose_level
+            if(rx[i]) if_debug_level
             {
-                double m = min(epsilon, beta[i]);
+                const double m = fmin(epsilon, beta[i]);
                 const double *r = &rays[nrows * i];
-                time_printf("  r[%d]=%.8lf %.8lf\n", i, r[0], r[1]);
+                time_printf("      r[%d]: %12.8lf %12.8lf\n", i, m * sbar[i] * r[0], m * sbar[i] * r[1]);
-                time_printf("  r[%d]=%.8lf %.8lf\n", i, m * r[0], m * r[1]);
             }
         }
     }
@@ -978,15 +1009,12 @@ static int bound(int nrows,
                 fbar[i] -= min(*epsilon, beta[j]) * r[i] * sbar[j];
         }
 
-        log_verbose("%.12lf %.12lf\n", f[0], f[1]);
-        log_verbose("%.12lf %.12lf\n", fbar[0], fbar[1]);
-
         prev_epsilon = *epsilon;
 
         rval = cone_bound(nrows, nrays, fbar, rays, rx, x, beta, epsilon);
         abort_if(rval, "cone_bound failed");
 
-        log_verbose("        e=%.12lf\n", *epsilon);
+        log_debug("    cone bound: %.20lf\n", epsilon);
         abort_if(prev_epsilon < *epsilon, "epsilon should never increase");
     }
 
@@ -1215,8 +1243,6 @@ CLEANUP:
 int INFINITY_ND_generate_lfree(const struct MultiRowModel *model,
                                struct ConvLFreeSet *lfree)
 {
-    log_debug("INFINITY_ND_generate_lfree\n");
-
     int rval = 0;
     int nrows = model->nrows;
     int nrays = model->rays.nrays;
@@ -1249,24 +1275,6 @@ int INFINITY_ND_generate_lfree(const struct MultiRowModel *model,
 
     int it = 0;
 
-    //lp_time = 0;
-    //lp_count = 0;
-
-    //epsilon_lp_count = 0;
-    //epsilon_lp_time = 0;
-    //
-    //sfree_mip_count = 0;
-    //sfree_mip_time = 0;
-
-    //tight_lp_count = 0;
-    //tight_lp_time = 0;
-
-    //violated_lp_count = 0;
-    //violated_lp_time = 0;
-
-    //scale_ahull_lp_time = 0;
-    //scale_ahull_lp_count = 0;
-
     long x_count = 0;
     double epsilon;
 
@@ -1289,7 +1297,7 @@ int INFINITY_ND_generate_lfree(const struct MultiRowModel *model,
 
         while(1)
         {
-            log_debug("  epsilon = %.8lf\n", epsilon);
+            log_debug("    epsilon: %.8lf\n", epsilon);
 
             int found = 0;
 
@@ -1309,8 +1317,14 @@ int INFINITY_ND_generate_lfree(const struct MultiRowModel *model,
                 if(!found) break;
             }
 
-            log_debug("    found interior point:\n");
+            if_debug_level {
-            for(int i = 0; i < nrows; i++) log_debug("        %.2lf\n", x[i]);
+                time_printf("    found interior point:\n");
+                time_printf("        x: ");
+                for (int i = 0; i < nrows; i++)
+                    printf("%.2lf ", x[i]);
+                printf("\n");
+
+            }
 
             x_count++;
             abort_if(x_count > 1000, "infinite loop");
@@ -1320,10 +1334,8 @@ int INFINITY_ND_generate_lfree(const struct MultiRowModel *model,
                     &epsilon_x, tx);
             abort_if(rval, "bound failed");
 
-            if(isinf(epsilon_x)) break;
+            abort_if(isinf(epsilon_x), "epsilon_x is infinite");
-
+            log_debug("    epsilon_x: %.8lf\n", epsilon_x);
-//            epsilon_x = (floor(epsilon_x * 128) / 128);
-            log_debug("    epsilon_x = %.8lf\n", epsilon_x);
 
             if(DOUBLE_eq(epsilon_x, epsilon))
             {
@@ -1343,11 +1355,12 @@ int INFINITY_ND_generate_lfree(const struct MultiRowModel *model,
 
         int skip_ahull = 0;
 
+        log_debug("    updating beta:\n");
         for(int i = 0; i < nrays; i++)
         {
             if(t[i])
             {
-                beta[i] = min(beta[i], epsilon);
+                beta[i] = min(beta[i] * 0.999, epsilon);
             }
             else if(!skip_ahull)
             {
@@ -1364,14 +1377,11 @@ int INFINITY_ND_generate_lfree(const struct MultiRowModel *model,
                     continue;
                 }
 
-//                alpha = (floor(alpha * 128) / 128);
                 beta[i] = min(beta[i], alpha);
             }
 
-            log_debug("  beta[%2d] = %20.12lf\n", i, beta[i]);
+            log_debug("        beta[%2d]: %20.12lf\n", i, beta[i]);
         }
-
-        log_debug("epsilon = %.6lf\n", epsilon);
     }
 
     log_debug("    %6ld lattice points, %ld iterations\n", x_count, it);

infinity/library/src/infinity.c

@@ -108,6 +108,7 @@ static int create_cut_from_lfree(const struct Tableau *tableau,
                                 struct Row *cut)
 {
     int rval = 0;
+    double *ray = 0;
 
     struct LP lp;
     int nvars = map->nvars;
@@ -122,40 +123,41 @@ static int create_cut_from_lfree(const struct Tableau *tableau,
     rval = INFINITY_create_psi_lp(lfree, &lp);
     abort_if(rval, "create_psi_lp failed");
 
+    ray = (double*) malloc(nrows * sizeof(double));
+    abort_if(!ray, "could not allocate ray");
+
     cut->nz = nvars;
     for(int i = 0; i < nvars; i++)
     {
-        double value;
+        double value, norm = 0;
-        const double *q = LFREE_get_ray(rays, map->variable_to_ray[i]);
+        double *original_ray = LFREE_get_ray(rays, map->variable_to_ray[i]);
         char type = tableau->column_types[map->indices[i]];
 
+        for (int j = 0; j < nrows; j++)
+        {
+            ray[j] = original_ray[j];
+            norm += fabs(ray[j]);
+        }
+
+        if (norm < 0.001)
+            for (int j = 0; j < nrows; j++)
+                ray[j] *= 0.001 / norm;
 
         if(ENABLE_LIFTING && type == MILP_INTEGER)
         {
-            rval = INFINITY_pi(nrows, q, map->ray_scale[i], &lp, &value);
+            rval = INFINITY_pi(nrows, ray, map->ray_scale[i], &lp, &value);
             abort_if(rval, "INFINITY_pi failed");
         }
         else
         {
-            rval = INFINITY_psi(nrows, q, map->ray_scale[i], &lp, &value);
+            rval = INFINITY_psi(nrows, ray, map->ray_scale[i], &lp, &value);
             abort_if(rval, "INFINITY_psi failed");
         }
 
+        value *= 1.001;
+        value = fmax(value, 0.001);
         log_verbose("   psi[%4d] = %20.12lf %d\n", map->indices[i], value);
 
-        if_debug_level
-        {
-            time_printf("    q=[");
-            for(int j = 0; j < nrows; j++)
-                printf("%25.20lf ", q[j] * map->ray_scale[i]);
-            printf("] value=%25.20lf\n", value);
-        }
-
-        value = fmax(value, 1 / 1024.0);
-
-//        value *= 1.001;
-//        value = DOUBLE_max(value, 0.001);
-
         cut->indices[i] = map->indices[i];
         cut->pi[i] = -value;
     }
@@ -163,6 +165,7 @@ static int create_cut_from_lfree(const struct Tableau *tableau,
     cut->pi_zero = -1.0;
 
 CLEANUP:
+    if(ray) free(ray);
     LP_free(&lp);
     return rval;
 }
@@ -186,12 +189,7 @@ static int filter_model(const struct MultiRowModel *original_model,
     struct RayList *filtered_rays = &filtered_model->rays;
     const struct RayList *original_rays = &original_model->rays;
 
-    memcpy(f, original_model->f, 2 * sizeof(double));
+    memcpy(f, original_model->f, nrows * sizeof(double));
-    for(int i = 0; i < nrows; i++)
-    {
-        f[i] = (ceil(f[i] * 128) / 128);
-        if(f[i] <= 0.01) f[i] = 0;
-    }
 
     r = (double*) malloc(nrows * sizeof(double));
     abort_if(!r, "could not allocate r");
@@ -268,11 +266,18 @@ static int append_extra_rays(struct MultiRowModel *model)
 
     for(int i = 0; i < nrows; i++)
     {
+        int found, index;
+        double scale;
+
         for(int j = 0; j < nrows; j++) r[j] = (i == j ? e : 0);
-        LFREE_push_ray(&model->rays, r);
+        rval = CG_find_ray(&model->rays, r, &found, &scale, &index);
+        abort_if(rval, "CG_find_ray failed");
+        if(!found) LFREE_push_ray(&model->rays, r);
 
         for(int j = 0; j < nrows; j++) r[j] = (i == j ? -e : 0);
-        LFREE_push_ray(&model->rays, r);
+        rval = CG_find_ray(&model->rays, r, &found, &scale, &index);
+        abort_if(rval, "CG_find_ray failed");
+        if(!found) LFREE_push_ray(&model->rays, r);
     }
 
     if(nrows == 2)
@@ -409,6 +414,9 @@ int INFINITY_generate_cut(const struct Tableau *tableau, struct Row *cut)
     rval = CG_init_map(&original_map, max_nrays, tableau->nrows);
     abort_if(rval, "CG_init_map failed");
 
+    rval = LFREE_init_conv(&lfree, tableau->nrows, max_nrays);
+    abort_if(rval, "LFREE_init_conv failed");
+
     rval = extract_models(tableau, &original_model, &filtered_model,
             &original_map);
     abort_if(rval, "extract_models failed");
@@ -419,41 +427,18 @@ int INFINITY_generate_cut(const struct Tableau *tableau, struct Row *cut)
         goto CLEANUP;
     }
 
-    rval = LFREE_init_conv(&lfree, tableau->nrows, max_nrays);
-    abort_if(rval, "LFREE_init_conv failed");
-
-//    if(tableau->nrows == 2)
-//        rval = INFINITY_2D_generate_lfree(&filtered_model, &lfree);
-//    else
-//        rval = INFINITY_ND_generate_lfree(&filtered_model, &lfree);
-//
-//    if(rval)
-//    {
-//        rval = ERR_NO_CUT;
-//        goto CLEANUP;
-//    }
-//
-//    for(int i = 0; i < filtered_model.rays.nrays; i++)
-//    {
-//        double *r = LFREE_get_ray(&filtered_model.rays, i);
-//        for(int j = 0; j < tableau->nrows; j++)
-//            r[j] *= lfree.beta[j];
-//    }
-
     rval = append_extra_rays(&filtered_model);
     abort_if(rval, "append_extra_rays failed");
 
-//    if(tableau->nrows == 2)
+    if_debug_level
-//        rval = INFINITY_2D_generate_lfree(&filtered_model, &lfree);
-//    else
-        rval = INFINITY_ND_generate_lfree(&filtered_model, &lfree);
-
-    if(rval)
     {
-        rval = ERR_NO_CUT;
+        rval = CG_print_model(&filtered_model);
-        goto CLEANUP;
+        abort_if(rval, "CG_print_model failed");
     }
 
+    rval = INFINITY_ND_generate_lfree(&filtered_model, &lfree);
+    abort_if(rval, "INFINITY_ND_generate_lfree failed");
+
     if(SHOULD_DUMP_CUTS)
     {
         rval = dump_cut(&lfree);

infinity/library/tests/infinity-nd-test.cpp

@@ -27,6 +27,7 @@ extern "C" {
 #include "../src/infinity-nd.c"
 }
 
+const double E = 1e-6;
 
 TEST(InfinityNDTest, find_violated_cone_test)
 {
@@ -164,11 +165,11 @@ TEST(InfinityNDTest, cone_bound_test_1)
 
     rval = cone_bound(2, 6, f, rays, rx1, x, beta, &epsilon);
     abort_if(rval, "cone_bound failed");
-    EXPECT_NEAR(0.5, epsilon, 1e-6);
+    EXPECT_NEAR(0.5, epsilon, E);
 
     rval = cone_bound(2, 6, f, rays, rx2, x, beta, &epsilon);
     abort_if(rval, "cone_bound failed");
-    EXPECT_NEAR(1.0, epsilon, 1e-6);
+    EXPECT_NEAR(1.0, epsilon, E);
 
     CLEANUP:
     if(rval) FAIL();
@@ -197,7 +198,7 @@ TEST(InfinityNDTest, cone_bound_test_2)
 
     rval = cone_bound(2, 2, f, rays, rx, x1, beta1, &epsilon);
     abort_if(rval, "cone_bound failed");
-    EXPECT_NEAR(1.0, epsilon, 1e-6);
+    EXPECT_NEAR(1.0, epsilon, E);
 
     rval = cone_bound(2, 2, f, rays, rx, x1, beta2, &epsilon);
     abort_if(rval, "cone_bound failed");
@@ -205,7 +206,7 @@ TEST(InfinityNDTest, cone_bound_test_2)
 
     rval = cone_bound(2, 2, f, rays, rx, x2, beta2, &epsilon);
     abort_if(rval, "cone_bound failed");
-    EXPECT_NEAR(1.0, epsilon, 1e-6);
+    EXPECT_NEAR(1.0, epsilon, E);
 
     rval = cone_bound(2, 2, f, rays, rx, x2, beta3, &epsilon);
     abort_if(rval, "cone_bound failed");
@@ -239,7 +240,7 @@ TEST(InfinityNDTest, bound_test_1)
 
     rval = bound(2, 6, f, rays, x, beta1, &epsilon, tx);
     abort_if(rval, "bound failed");
-    EXPECT_NEAR(epsilon, 0.5, 1e-6);
+    EXPECT_NEAR(epsilon, 0.5, E);
     EXPECT_TRUE(tx[0]);
     EXPECT_FALSE(tx[1]);
     EXPECT_FALSE(tx[2]);
@@ -249,7 +250,7 @@ TEST(InfinityNDTest, bound_test_1)
 
     rval = bound(2, 6, f, rays, x, beta2, &epsilon, tx);
     abort_if(rval, "bound failed");
-    EXPECT_NEAR(epsilon, 1.0, 1e-6);
+    EXPECT_NEAR(epsilon, 1.0, E);
     EXPECT_TRUE(tx[0]);
     EXPECT_FALSE(tx[1]);
     EXPECT_FALSE(tx[2]);
@@ -308,11 +309,11 @@ TEST(InfinityNDTest, psi_test)
 
     rval = INFINITY_psi(2, q1, 1.0, &lp, &value);
     abort_if(rval, "GREDDY_ND_psi failed");
-    EXPECT_NEAR(value, 2.0, 1e-6);
+    EXPECT_NEAR(value, 2.0, E);
 
     rval = INFINITY_psi(2, q2, 2.0, &lp, &value);
     abort_if(rval, "GREDDY_ND_psi failed");
-    EXPECT_NEAR(value, 8.0, 1e-6);
+    EXPECT_NEAR(value, 8.0, E);
 
 CLEANUP:
     LP_free(&lp);
@@ -360,17 +361,64 @@ TEST(InfinityNDTest, psi_test_2)
 
     rval = INFINITY_psi(3, q1, 1.0, &lp, &value);
     abort_if(rval, "GREDDY_ND_psi failed");
-    EXPECT_NEAR(value, 1.0, 1e-6);
+    EXPECT_NEAR(value, 1.0, E);
 
     rval = INFINITY_psi(3, q2, 1.0, &lp, &value);
     abort_if(rval, "GREDDY_ND_psi failed");
-    EXPECT_NEAR(value, 2.0, 1e-6);
+    EXPECT_NEAR(value, 2.0, E);
 
 CLEANUP:
     if(rval) FAIL();
 }
 
-TEST(InfinityNDTest, generate_cut_test_1)
+TEST(InfinityNDTest, psi_test_3)
+{
+    int rval = 0;
+
+    double f[] = { 0.671875, 0.671875 };
+    double rays[] = {
+            -0.007812500000,  0.000000000000,
+            -0.039062500000,  0.046875000000,
+             0.000000000000,  0.046875000000,
+             0.046875000000,  0.000000000000,
+             0.000000000000, -0.039062500000
+    };
+
+    double beta[] = {
+            66.909090909091,
+            29.440000000000,
+            14.000000000000,
+            14.000000000000,
+            29.440000000000,
+    };
+
+    double q[] = { 0 - f[0], 1 - f[1]};
+
+    struct ConvLFreeSet lfree;
+    lfree.f = f;
+    lfree.beta = beta;
+    lfree.rays.nrays = 5;
+    lfree.rays.values = rays;
+    lfree.nrows = lfree.rays.dim = 2;
+
+    double value;
+    struct LP lp;
+
+    rval = LP_open(&lp);
+    abort_if(rval, "LP_open failed");
+
+    rval = INFINITY_create_psi_lp(&lfree, &lp);
+    abort_if(rval, "INFINITY_create_psi_lp failed");
+
+    rval = INFINITY_psi(lfree.nrows, q, 1.0, &lp, &value);
+    abort_if(rval, "GREDDY_ND_psi failed");
+    EXPECT_NEAR(value, 1.0, E);
+
+CLEANUP:
+    if(rval) FAIL();
+}
+
+TEST(DISABLED_InfinityNDTest, generate_cut_test_1)
 {
     int rval = 0;
 
@@ -398,12 +446,12 @@ TEST(InfinityNDTest, generate_cut_test_1)
     rval = INFINITY_ND_generate_lfree(&model, &lfree);
     abort_if(rval, "INFINITY_ND_generate_lfree failed");
 
-    EXPECT_NEAR(lfree.beta[0], 0.5, 1e-6);
+    EXPECT_NEAR(lfree.beta[0], 0.5, E);
-    EXPECT_NEAR(lfree.beta[1], 0.5, 1e-6);
+    EXPECT_NEAR(lfree.beta[1], 0.5, E);
-    EXPECT_NEAR(lfree.beta[2], 0.5, 1e-6);
+    EXPECT_NEAR(lfree.beta[2], 0.5, E);
-    EXPECT_NEAR(lfree.beta[3], 0.5, 1e-6);
+    EXPECT_NEAR(lfree.beta[3], 0.5, E);
-    EXPECT_NEAR(lfree.beta[4], 1.0, 1e-6);
+    EXPECT_NEAR(lfree.beta[4], 1.0, E);
-    EXPECT_NEAR(lfree.beta[5], 1.0, 1e-6);
+    EXPECT_NEAR(lfree.beta[5], 1.0, E);
 
 CLEANUP:
     if(rval) FAIL();
@@ -438,12 +486,12 @@ TEST(InfinityNDTest, generate_cut_test_2)
     rval = INFINITY_ND_generate_lfree(&model, &lfree);
     abort_if(rval, "INFINITY_ND_generate_lfree failed");
 
-    EXPECT_NEAR(lfree.beta[0], 0.75, 1e-6);
+    EXPECT_NEAR(lfree.beta[0], 0.75, E);
-    EXPECT_NEAR(lfree.beta[1], 2.25, 1e-6);
+    EXPECT_NEAR(lfree.beta[1], 2.25, E);
-    EXPECT_NEAR(lfree.beta[2], 0.75, 1e-6);
+    EXPECT_NEAR(lfree.beta[2], 0.75, E);
-    EXPECT_NEAR(lfree.beta[3], 2.25, 1e-6);
+    EXPECT_NEAR(lfree.beta[3], 2.25, E);
-    EXPECT_NEAR(lfree.beta[4], 0.75, 1e-6);
+    EXPECT_NEAR(lfree.beta[4], 0.75, E);
-    EXPECT_NEAR(lfree.beta[5], 2.25, 1e-6);
+    EXPECT_NEAR(lfree.beta[5], 2.25, E);
 
 CLEANUP:
     CG_free_model(&model);
@@ -487,3 +535,68 @@ TEST(InfinityNDTest, scale_to_ahull_test)
 CLEANUP:
     if(rval) FAIL();
 }
+
+TEST(InfinityNDTest, cone_bound_find_lambda_test)
+{
+    int rval = 0;
+
+    double f[] = { 1/2.0, 1/3.0 };
+    double x[] = { 3.0, 3.0 };
+    double ray_values[] = {
+            1.0, 2.0,
+            3.0, 2.0,
+    };
+
+    struct RayList rays = {.values = ray_values, .nrays = 2, .dim = 2};
+    double lambda[2];
+
+    rval = cone_bound_find_lambda(&rays, f, x, lambda);
+    abort_if(rval, "cone_bound_find_lambda failed");
+
+    EXPECT_NEAR(lambda[0], 3/4.0, E);
+    EXPECT_NEAR(lambda[1], 7/12.0, E);
+
+ CLEANUP:
+    if(rval) FAIL();
+}
+
+TEST(InfinityNDTest, cone_bound_find_lambda_test_2)
+{
+    int rval = 0;
+
+    double f[] = { 1/2.0, 1/3.0, 5/6.0 };
+    double x[] = { 3.0, 3.0, 6.0 };
+    double ray_values[] = {
+            1.0, 2.0, 3.0,
+            3.0, 2.0, 5.0,
+    };
+
+    struct RayList rays = {.values = ray_values, .nrays = 2, .dim = 3};
+    double lambda[2];
+
+    rval = cone_bound_find_lambda(&rays, f, x, lambda);
+    abort_if(rval, "cone_bound_find_lambda failed");
+
+    EXPECT_NEAR(lambda[0], 3/4.0, E);
+    EXPECT_NEAR(lambda[1], 7/12.0, E);
+
+CLEANUP:
+    if(rval) FAIL();
+}
+
+TEST(InfinityNDTest, cone_bound_find_lambda_test_3)
+{
+    int rval = 0;
+    double f[] = { 1/2.0, 1/2.0 };
+    double x[] = { 3.0, 3.0 };
+    double ray_values[] = {
+             0.0, -1.0,
+            -1.0,  0.0,
+    };
+
+    struct RayList rays = {.values = ray_values, .nrays = 2, .dim = 2};
+    double lambda[2];
+
+    rval = cone_bound_find_lambda(&rays, f, x, lambda);
+    EXPECT_NE(rval, 0);
+}

lifting/README.md

@@ -0,0 +1,52 @@
+lifting
+=======
+
+This package contains the source code for the paper **The (not so) Trivial
+Lifting in Two Dimensions**, by Ricardo Fukasawa, Laurent Poirrier and Álinson
+S. Xavier.
+
+Required Tools and Libraries
+----------------------------
+
+To produce the tables in the paper, the following tools and libraries were
+used. Different versions may produce slightly different outputs.
+
+- GNU Make, version 3.81
+- CMake, version 3.7.2
+- GCC, the GNU Compiler Collection, version 6.3.0
+- Ruby, version 2.4.0
+- IBM® ILOG® CPLEX®, version 12.6
+
+Build instructions
+------------------
+
+1. Navigate to the folder `../build`
+2. Run `cmake ..` followed by `make lifting-benchmark.run`
+3. Two binaries (`lifting-benchmark.run` and `liblifting.a`) will be generated
+
+Running the experiments
+-----------------------
+
+1. Build the project, following the instructions above.
+2. Navigate to the folder `lifting/benchmark` and execute
+
+        ./run_experiments.sh
+
+3. Two CSV files will be generated inside the folder `lifting/benchmark/tables`,
+	corresponding to the two tables that appear in the paper.
+
+Modifying the instances
+-----------------------
+
+In order to run the experiments with a different set of instances,
+the file `lifting/benchmark/instances/filtered/all.txt` should be modified.
+Each line in this file describes the origin `f` and a lattice-free set `B`.
+The set `B` is described by the coordinates of its vertices. Since the
+benchmark code only deals with maximal lattice-free sets, it is also necessary
+to specify the lattice-points that belong to each facet of `B`.
+If `n` is the number of facets, `v` is an n-by-2 matrix of doubles corresponding
+to the vertices and l is an n-by-2 matrix of doubles corresponding to the
+lattice-points, then each line of the file should be written as
+
+	f[0] f[1] n v[0][0] v[0][1] ... v[n-1][0] v[n-1][1] n l[0][0] l[0][1] ... l[n-1][0] l[n-1][1]
+

lifting/benchmark/answers/shear-1240.txt

@@ -6917,7 +6917,7 @@
 69 16 0.68683156580163995386
 69 17 0.92943801855653873645
 69 18 0.10740586958127096295
-69 19 0.92496934696100652218
+69 19 0.92497488891240209341
 69 20 0.60183083702384010394
 69 21 0.49150980007834732533
 69 22 0.42780977682599541367
@@ -51695,7 +51695,7 @@
 516 94 0.76878394367395230802
 516 95 0.38042760308826473192
 516 96 0.39931448388506396441
-516 97 0.47279282059935212601
+516 97 0.47279361768801209109
 516 98 0.97282870978187929722
 516 99 0.62348135864249343285
 517 0 0.50616530635928924653
@@ -94349,7 +94349,7 @@
 943 48 0.76196620377269397295
 943 49 0.22952707635243996265
 943 50 0.41734143911918408776
-943 51 0.91839681508463399950
+943 51 0.91840633855952091835
 943 52 0.37241285013528369063
 943 53 0.82926766528856887817
 943 54 0.54408107882613876427

lifting/benchmark/run_experiments.sh

@@ -8,6 +8,8 @@ function title()
 }
 
 RUN=../../build/lifting/benchmark/lifting-benchmark.run
+make -C ../../build lifting-benchmark.run || exit 1
+
 if [ ! -f $RUN ]; then
     echo "not found: $RUN"
     echo "please build the project before running this script"
@@ -15,12 +17,6 @@ if [ ! -f $RUN ]; then
 fi
 
 INSTANCES="instances/filtered/all.txt"
-# SAMPLES_SLOW=10
-# SAMPLES_MEDIUM=100
-# SAMPLES_FAST=1000
-SAMPLES_SLOW=1
-SAMPLES_MEDIUM=1
-SAMPLES_FAST=1
 SEED=1240
 
 # ORIGINAL
@@ -36,21 +32,30 @@ COMMON_OPTS="$COMMON_OPTS --check-answers $ANSWERS"
 DIR=orig-100
 mkdir -p $DIR; rm -f $DIR/*log $DIR/*yaml
 
+title Heuristic
+$RUN $COMMON_OPTS --samples 100000 --heuristic --log $DIR/heur.log --stats $DIR/heur.yaml || exit
+
 title Bound Original
-$RUN $COMMON_OPTS --samples $SAMPLES_FAST --bound --log $DIR/bound-nopre.log --stats $DIR/bound-nopre.yaml || exit
+$RUN $COMMON_OPTS --samples 1000 --bound --log $DIR/bound.log --stats $DIR/bound.yaml || exit
 
 title Bound Pre-processing
-$RUN $COMMON_OPTS --samples $SAMPLES_FAST --bound --preprocess --log $DIR/bound-pre.log --stats $DIR/bound-pre.yaml || exit
+$RUN $COMMON_OPTS --samples 1000 --bound --preprocess --log $DIR/bound-pre.log --stats $DIR/bound-pre.yaml || exit
 
 title Naive Bounding-Box
-$RUN $COMMON_OPTS --samples $SAMPLES_MEDIUM --naive --log $DIR/naive-bbox.log --stats $DIR/naive-bbox.yaml || exit
+$RUN $COMMON_OPTS --samples 100 --naive --log $DIR/naive-bbox.log --stats $DIR/naive-bbox.yaml || exit
 
 title Naive Fixed-M
 M=50
-$RUN $COMMON_OPTS --samples $SAMPLES_MEDIUM --naive --fixed-bounds $M --log $DIR/naive-fixed-$M.log --stats $DIR/naive-fixed-$M.yaml || exit
+$RUN $COMMON_OPTS --samples 100 --naive --fixed-bounds $M --log $DIR/naive-fixed-$M.log --stats $DIR/naive-fixed-$M.yaml || exit
+
+title Naive Bounding-Box Pre-processing
+$RUN $COMMON_OPTS --samples 100 --naive --preprocess --log $DIR/naive-bbox-pre.log --stats $DIR/naive-bbox-pre.yaml || exit
 
 title MIP
-$RUN $COMMON_OPTS --samples $SAMPLES_SLOW --mip --log $DIR/mip.log --stats $DIR/mip.yaml || exit
+$RUN $COMMON_OPTS --samples 10 --mip --log $DIR/mip.log --stats $DIR/mip.yaml || exit
+
+title MIP Pre-processing
+$RUN $COMMON_OPTS --samples 10 --mip --preprocess --log $DIR/mip-pre.log --stats $DIR/mip-pre.yaml || exit
 
 # SHEAR
 # ------------------------------------------------------------------------------
@@ -65,21 +70,28 @@ COMMON_OPTS="$COMMON_OPTS --check-answers $ANSWERS"
 DIR=shear-100
 mkdir -p $DIR; rm -f $DIR/*log $DIR/*yaml
 
-title Bound Pre-processing + Shear
+title Heuristic + Shear
-$RUN $COMMON_OPTS --samples $SAMPLES_FAST --bound --preprocess --log $DIR/bound-pre.log --stats $DIR/bound-pre.yaml || exit
+$RUN $COMMON_OPTS --samples 100000 --heuristic --log $DIR/heur.log --stats $DIR/heur.yaml || exit
 
 title Bound Original + Shear
-$RUN $COMMON_OPTS --samples $SAMPLES_MEDIUM --bound --log $DIR/bound-nopre.log --stats $DIR/bound-nopre.yaml || exit
+$RUN $COMMON_OPTS --samples 1000 --bound --log $DIR/bound.log --stats $DIR/bound.yaml || exit
+
+title Bound Pre-processing + Shear
+$RUN $COMMON_OPTS --samples 1000 --bound --preprocess --log $DIR/bound-pre.log --stats $DIR/bound-pre.yaml || exit
 
 title Naive Bounding-Box + Shear
-$RUN $COMMON_OPTS --samples $SAMPLES_SLOW --naive --log $DIR/naive-bbox.log --stats $DIR/naive-bbox.yaml || exit
+$RUN $COMMON_OPTS --samples 10 --naive --log $DIR/naive-bbox.log --stats $DIR/naive-bbox.yaml || exit
 
 title Naive Fixed-M + Shear
 M=50
-$RUN $COMMON_OPTS --samples $SAMPLES_MEDIUM --naive --fixed-bounds $M --log $DIR/naive-fixed-$M.log --stats $DIR/naive-fixed-$M.yaml || exit
+$RUN $COMMON_OPTS --samples 100 --naive --fixed-bounds $M --log $DIR/naive-fixed-$M.log --stats $DIR/naive-fixed-$M.yaml || exit
+
+title MIP + Shear
+$RUN $COMMON_OPTS --samples 10 --mip --log $DIR/mip.log --stats $DIR/mip.yaml || exit
+
+title MIP Pre-processing + Shear
+$RUN $COMMON_OPTS --samples 10 --mip --preprocess --log $DIR/mip-pre.log --stats $DIR/mip-pre.yaml || exit
 
-title MIP
-$RUN $COMMON_OPTS --samples $SAMPLES_SLOW --mip --log $DIR/mip.log --stats $DIR/mip.yaml || exit
 
 # TABLES
 # ------------------------------------------------------------------------------

lifting/benchmark/scripts/table.rb

@@ -3,8 +3,6 @@
 require 'yaml'
 require 'bigdecimal'
 
-BIG_M = 1000000
-
 def sum(a)
     a.inject(0){ |accum, i| accum + i }
 end
@@ -35,7 +33,7 @@ ARGV.each_with_index do |filename,idx|
     filenames[idx] = filename
             .gsub(/.*\//,"")
             .gsub(".yaml", "")
-            .gsub(/[^A-Za-z]/, "")
+            .gsub(/[^A-Za-z0-9]/, "")
     files[idx] = YAML::load(File.open(filename))
 end
 
@@ -59,14 +57,14 @@ files.each_with_index do |f,idx|
     for i in 0..(n_instances-1) do
         time = f['cpu_time'][i]
         next if time.nil?
-
+        times[idx].push(time)
         all_times[idx].push(time)
 
-        if(time == BIG_M)
+        wrong = f['wrong_answer'][i]
+        if(wrong == 1)
             fail_count[idx] += 1
         else
             success_count[idx] += 1
-            times[idx].push(time)
         end
     end
 
@@ -75,6 +73,7 @@ end
 
 best_percentage = [0] * files.length
 
+BIG_M = 1000000000
 for i in 0..(n_instances-1) do
     best_time = BIG_M
     files.each_with_index do |f,idx|
@@ -87,7 +86,6 @@ for i in 0..(n_instances-1) do
         best_percentage[idx] += 1 if time == best_time
 
         next if best_time <= 0
-        next if time >= BIG_M
         ratios_to_best[idx].push(time / best_time)
     end
 end

lifting/benchmark/src/main.c

@@ -40,10 +40,12 @@ unsigned int SEED = 0;
 #define ALGORITHM_BOUND 0
 #define ALGORITHM_NAIVE 1
 #define ALGORITHM_MIP 2
+#define ALGORITHM_HEUR 3
 
 int SELECT_NAIVE_ALGORITHM = 0;
 int SELECT_BOUND_ALGORITHM = 0;
 int SELECT_MIP_ALGORITHM = 0;
+int SELECT_HEUR_ALGORITHM = 0;
 
 int ENABLE_PREPROCESSING = 0;
 int ENABLE_SHEAR = 0;
@@ -83,6 +85,7 @@ static const struct option options_tab[] =
     {"check-answers", required_argument, 0, 'c'},
     {"samples", required_argument, 0, 'a'},
     {"mip", no_argument, 0, 'm'},
+    {"heuristic", no_argument, 0, 'r'},
     {0, 0, 0, 0}
 };
 
@@ -108,6 +111,8 @@ static void print_usage(char **argv)
             "select MIP algorithm");
     printf("%4s %-20s %s\n", "-u", "--bound",
             "select bound algorithm");
+    printf("%4s %-20s %s\n", "-e", "--heuristic",
+            "select heuristic algorithm");
     printf("%4s %-20s %s\n", "-p", "--preprocess",
             "enable pre-processing step in bound algorithm");
     printf("%4s %-20s %s\n", "-e", "--shear",
@@ -144,7 +149,7 @@ static int parse_args(int argc,
     {
         int c = 0;
         int option_index = 0;
-        c = getopt_long(argc, argv, "hb:k:s:f:o:nupew:c:a:m", options_tab,
+        c = getopt_long(argc, argv, "hb:k:s:f:o:nupew:c:a:mr", options_tab,
                         &option_index);
 
         if (c < 0) break;
@@ -208,6 +213,10 @@ static int parse_args(int argc,
             SELECT_BOUND_ALGORITHM = 1;
             break;
 
+        case 'r':
+            SELECT_HEUR_ALGORITHM = 1;
+            break;
+
         case 'p':
             ENABLE_PREPROCESSING = 1;
             break;
@@ -243,7 +252,8 @@ static int parse_args(int argc,
         rval = 1;
     }
 
-    if (SELECT_NAIVE_ALGORITHM + SELECT_BOUND_ALGORITHM + SELECT_MIP_ALGORITHM != 1)
+    if (SELECT_NAIVE_ALGORITHM + SELECT_BOUND_ALGORITHM + SELECT_MIP_ALGORITHM
+            + SELECT_HEUR_ALGORITHM != 1)
     {
         fprintf(stderr, "You must select exactly one algorithm.\n");
         rval = 1;
@@ -298,6 +308,18 @@ int benchmark_set_sample(int algorithm,
                          int *wrong_answer)
 {
     int rval = 0;
+    double xi_plus, xi_minus, ignored;
+
+    if(algorithm == ALGORITHM_BOUND)
+    {
+        rval = LIFTING_2D_optimize_continuous(set->n_halfspaces,
+                set->halfspaces, 1, &ignored, &xi_plus);
+        abort_if(rval, "LIFTING_2D_optimize_continuous failed");
+
+        rval = LIFTING_2D_optimize_continuous(set->n_halfspaces,
+                set->halfspaces, -1, &ignored, &xi_minus);
+        abort_if(rval, "LIFTING_2D_optimize_continuous failed");
+    }
 
     for (int i = 0; i < N_RAYS; i++)
     {
@@ -319,8 +341,8 @@ int benchmark_set_sample(int algorithm,
         switch (algorithm)
         {
             case ALGORITHM_BOUND:
-                rval = LIFTING_2D_bound(set->n_halfspaces, set->halfspaces, ray,
+                rval = LIFTING_2D_bound(set->n_halfspaces, set->halfspaces,
-                        &value);
+                        ray, xi_plus, xi_minus, &value);
                 abort_if(rval, "LIFTING_2D_bound failed");
                 break;
 
@@ -336,6 +358,12 @@ int benchmark_set_sample(int algorithm,
                 abort_if(rval, "LIFTING_2D_mip failed");
                 break;
 
+            case ALGORITHM_HEUR:
+                rval = LIFTING_2D_heur(set->n_halfspaces, set->halfspaces, ray,
+                        &value);
+                abort_if(rval, "LIFTING_2D_heur failed");
+                break;
+
             default:
                 abort_if(1, "Invalid algorithm");
         }
@@ -371,8 +399,6 @@ int benchmark_set_sample(int algorithm,
                             " expected=%.8lf delta=%.8lf)\n", set_idx,
                             i, value, expected_value, delta);
                     *wrong_answer = 1;
-
-                    LFREE_2D_print_set(set);
                 }
             }
 
@@ -440,11 +466,17 @@ int benchmark(int n_sets, struct LFreeSet2D *sets, double *rays,
         int algorithm)
 {
     int rval = 0;
+    double *times = 0;
+    int *wrong = 0;
+
+    wrong = (int*) malloc(n_sets * sizeof(int));
+    times = (double*) malloc(n_sets * sizeof(double));
+    abort_if(!wrong, "could not allocate wrong");
+    abort_if(!times, "could not allocate times");
 
     log_info("Running benchmark...\n");
 
     double total_initial_time = get_user_time();
-    stats_printf("cpu_time:\n");
 
     for (int j = 0; j < n_sets; j++)
     {
@@ -467,12 +499,20 @@ int benchmark(int n_sets, struct LFreeSet2D *sets, double *rays,
         double set_duration = get_user_time() - set_initial_time;
         double avg = (set_duration / N_SAMPLES_PER_SET) * 1000;
 
-        if(wrong_answer) avg = 1000000;
+        times[j] = avg;
+        wrong[j] = wrong_answer;
 
-        stats_printf("  %d: %.8lf\n", j, avg);
         log_info("  %3d: %12.3lf ms\n", j, avg);
     }
 
+    stats_printf("cpu_time:\n");
+    for(int j = 0; j < n_sets; j++)
+        stats_printf("  %d: %.8lf\n", j, times[j]);
+
+    stats_printf("wrong_answer:\n");
+    for(int j = 0; j < n_sets; j++)
+        stats_printf("  %d: %d\n", j, wrong[j]);
+
     double total_duration = get_user_time() - total_initial_time;
 
     log_info("    %.3lf ms per set                     \n",
@@ -485,6 +525,8 @@ int benchmark(int n_sets, struct LFreeSet2D *sets, double *rays,
     }
 
 CLEANUP:
+    if(wrong) free(wrong);
+    if(times) free(times);
     return rval;
 }
 
@@ -550,22 +592,32 @@ int main(int argc, char **argv)
         ray[1] = DOUBLE_random(0.0, 1.0);
     }
 
-    int algorithm = ALGORITHM_BOUND;
+    int algorithm = -1;
-    if(SELECT_MIP_ALGORITHM) algorithm = ALGORITHM_MIP;
-    else if(SELECT_NAIVE_ALGORITHM) algorithm = ALGORITHM_NAIVE;
     
-    if(algorithm == ALGORITHM_NAIVE)
+    if(SELECT_BOUND_ALGORITHM)
+    {
+        log_info("Enabling bound algorithm\n");
+        algorithm = ALGORITHM_BOUND;
+    }
+    else if(SELECT_MIP_ALGORITHM)
+    {
+        log_info("Enabling MIP algorithm\n");
+        algorithm = ALGORITHM_MIP;
+    }
+    else if(SELECT_NAIVE_ALGORITHM)
     {
         log_info("Enabling naive algorithm\n");
+        algorithm = ALGORITHM_NAIVE;
 
         if(USE_FIXED_BOUNDS)
             log_info("Using fixed big M: %d\n", NAIVE_BIG_M);
         else
             log_info("Enabling bounding boxes\n");
     }
-    else
+    else if(SELECT_HEUR_ALGORITHM)
     {
-        log_info("Enabling bound algorithm\n");
+        log_info("Enabling heuristic algorithm\n");
+        algorithm = ALGORITHM_HEUR;
     }
 
     log_info("Setting %d samples per set\n", N_SAMPLES_PER_SET);

lifting/benchmark/tables/orig-100.csv

@@ -1,7 +1,7 @@
-metric,boundnopre,boundpre,mip,naivebbox,naivefixed,
+metric,boundpre,bound,heur,mippre,mip,naivebboxpre,naivebbox,naivefixed50,
-Average (ms),0.565,0.059,141.938,9.341,23.278,
+Average (ms),0.078,0.937,0.004,70.080,81.738,6.353,12.490,34.874,
-Median (ms),0.068,0.060,129.200,2.120,23.080,
+Median (ms),0.079,0.099,0.004,67.590,73.889,2.900,3.120,34.195,
-Maximum (ms),33.264,0.100,947.600,1769.880,30.080,
+Maximum (ms),0.152,55.998,0.008,175.773,1705.441,96.585,2483.652,122.881,
-Failure Rate,0.0 \%,0.0 \%,0.0 \%,0.0 \%,0.3 \%,
+Failure Rate,0.0 \%,0.0 \%,100.0 \%,0.0 \%,0.0 \%,0.0 \%,0.1 \%,0.3 \%,
-Best,37.2 \%,83.6 \%,0.0 \%,0.0 \%,0.0 \%,
+Best,0.0 \%,0.0 \%,100.0 \%,0.0 \%,0.0 \%,0.0 \%,0.0 \%,0.0 \%,
-Avg Ratio to Best,9.483,1.015,2449.938,159.449,402.950,
+Avg Ratio to Best,18.918,223.183,1.000,17133.363,20066.439,1539.254,3223.466,8511.913,

lifting/benchmark/tables/shear-100.csv

@@ -1,7 +1,7 @@
-metric,boundnopre,boundpre,mip,naivebbox,naivefixed,
+metric,boundpre,bound,mip,naivebbox,naivefixed50,
-Average (ms),4.709,0.062,1906.081,5473.688,23.374,
+Average (ms),0.080,7.723,659.486,478.703,34.680,
-Median (ms),0.440,0.064,582.800,21.200,23.160,
+Median (ms),0.080,0.730,131.680,31.395,34.105,
-Maximum (ms),485.760,0.108,94620.000,392259.600,32.120,
+Maximum (ms),0.352,800.768,53099.128,9450.563,48.423,
-Failure Rate,0.0 \%,0.0 \%,1.7 \%,0.1 \%,10.0 \%,
+Failure Rate,0.0 \%,0.0 \%,2.1 \%,0.5 \%,10.0 \%,
-Best,0.0 \%,100.0 \%,0.0 \%,0.0 \%,0.0 \%,
+Best,100.0 \%,0.0 \%,0.0 \%,0.0 \%,0.0 \%,
-Avg Ratio to Best,74.602,1.000,31250.644,87472.010,380.532,
+Avg Ratio to Best,1.000,96.571,8243.616,6018.030,446.701,

lifting/library/include/lifting/lifting.h

@@ -46,8 +46,15 @@ int LIFTING_2D_naive(int n_halfspaces,
 int LIFTING_2D_bound(int n_halfspaces,
                      const double *halfspaces,
                      const double *ray,
+                     const double xi_plus,
+                     const double xi_minus,
                      double *value);
 
+int LIFTING_2D_heur(int n_halfspaces,
+                    const double *halfspaces,
+                    const double *ray,
+                    double *value);
+
 int LIFTING_2D_verify(struct LFreeSet2D *set);
 
 #endif //LIFTING_H

lifting/library/src/lifting.c

@@ -205,15 +205,32 @@ CLEANUP:
     return rval;
 }
 
+int LIFTING_2D_heur(int n_halfspaces,
+                    const double *halfspaces,
+                    const double *ray,
+                    double *value)
+{
+    int rval = 0;
+
+    double q[2] = { ray[0] - ceil(ray[0]), ray[1] - ceil(ray[1])};
+
+    rval = LIFTING_2D_psi(n_halfspaces, halfspaces, q, value);
+    abort_if(rval, "LIFTING_2D_ps failed");
+
+CLEANUP:
+    return rval;
+}
+
 int LIFTING_2D_bound(int n_halfspaces,
                      const double *halfspaces,
                      const double *ray,
+                     const double xi_plus,
+                     const double xi_minus,
                      double *value)
 {
     int rval = 0;
 
     double eta_star, eta_plus, eta_minus;
-    double xi_plus, xi_minus;
     double m_plus, m_minus;
     double ignored;
     int k1 = 1;
@@ -221,14 +238,6 @@ int LIFTING_2D_bound(int n_halfspaces,
     rval = LIFTING_2D_lift_fixed(n_halfspaces, halfspaces, ray, 0, &eta_star);
     abort_if(rval, "LIFTING_2D_lift_fixed failed");
 
-    rval = LIFTING_2D_optimize_continuous(n_halfspaces, halfspaces, 1, &ignored,
-            &xi_plus);
-    abort_if(rval, "LIFTING_2D_optimize_continuous failed");
-
-    rval = LIFTING_2D_optimize_continuous(n_halfspaces, halfspaces, -1,
-            &ignored, &xi_minus);
-    abort_if(rval, "LIFTING_2D_optimize_continuous failed");
-
     m_plus = m_minus = INFINITY;
 
     log_debug("Level 0:\n");

multirow/CMakeLists.txt

@@ -1,33 +1,43 @@
 set(COMMON_SOURCES
-    src/cg.c
+        src/cg.c
-    src/double.c
+        src/double.c
-    src/geometry.c
+        src/geometry.c
-    src/lfree2d.c
+        src/lfree2d.c
-    src/lp.c
+        src/lp.c
-    src/mir.c
+        src/mir.c
-    src/util.c
+        src/util.c
-    src/rational.c
+        src/rational.c
-    src/stats.c
+        src/stats.c
-    include/multirow/cg.h
+        src/linalg.c
-    include/multirow/double.h
+        include/multirow/cg.h
-    include/multirow/geometry.h
+        include/multirow/double.h
-    include/multirow/lfree2d.h
+        include/multirow/geometry.h
-    include/multirow/lp.h
+        include/multirow/lfree2d.h
-    include/multirow/mir.h
+        include/multirow/lp.h
-    include/multirow/rational.h
+        include/multirow/mir.h
-    include/multirow/stats.h
+        include/multirow/rational.h
-    include/multirow/params.h
+        include/multirow/stats.h
-    include/multirow/util.h)
+        include/multirow/params.h
+        include/multirow/util.h
+        include/multirow/linalg.h)
 
 set(TEST_SOURCES
-    tests/double-test.cpp
+        tests/double-test.cpp
-    tests/cg-test.cpp
+        tests/cg-test.cpp
-    tests/geometry-test.cpp)
+        tests/geometry-test.cpp
+        tests/linalg-test.cpp)
 
 add_library(multirow_static ${COMMON_SOURCES})
 set_target_properties(multirow_static PROPERTIES OUTPUT_NAME lifting)
-target_link_libraries(multirow_static ${CPLEX_LIBRARIES})
+target_link_libraries(multirow_static
-target_include_directories (multirow_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
+        ${CPLEX_LIBRARIES}
+        ${BLAS_LIBRARIES}
+        ${LAPACKE_LIBRARIES})
+target_include_directories(multirow_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 
 add_executable(multirow-test.run ${COMMON_SOURCES} ${TEST_SOURCES})
-target_link_libraries(multirow-test.run ${CPLEX_LIBRARIES} gtest_main)
+target_link_libraries(multirow-test.run
+        ${CPLEX_LIBRARIES}
+        ${BLAS_LIBRARIES}
+        ${LAPACKE_LIBRARIES}
+        gtest_main)

multirow/include/multirow/cg.h

@@ -116,4 +116,6 @@ int CG_total_nz(const struct Tableau *tableau);
 
 double CG_replace_x(const struct Row *row, const double *x);
 
+int CG_print_model(const struct MultiRowModel *model);
+
 #endif //MULTIROW_CG_H

multirow/include/multirow/lfree2d.h

@@ -40,7 +40,6 @@ struct RayList
     int dim;
 };
 
-
 struct ConvLFreeSet
 {
     double *f;
@@ -82,4 +81,6 @@ int LFREE_init_conv(struct ConvLFreeSet *lfree, int dim, int max_nrays);
 
 void LFREE_free_conv(struct ConvLFreeSet *lfree);
 
+int LFREE_print_set(const struct ConvLFreeSet *lfree);
+
 #endif //LFREE_2D_H

multirow/include/multirow/linalg.h

@@ -0,0 +1,40 @@
+/* Copyright (c) 2015-2017 Alinson Xavier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+/**
+ * Receives an n-dimensional vector x, a scalar alpha and sets x <- alpha * x
+ */
+void LINALG_scale(int n, double *x, double alpha);
+
+/**
+ * Receives two n-dimensional vectors x and y and returns the dot product of x
+ * and y.
+ */
+double LINALG_dot(int n, double *x, double *y);
+
+/**
+ * Receives an n-dimensional vector x and returns the 1-norm of x.
+ */
+double LINALG_norm(int n, double *x);
+
+/**
+ * Given a full rank m-by-n matrix A and an m-dimensional vector b, this
+ * function finds x such that Ax = b. Returns zero if the operation is
+ * successful and non-zero otherwise.
+ */
+int LINALG_solve(int n, int m, double *A, double *b, double *x);

multirow/include/multirow/params.h

@@ -36,7 +36,7 @@
 /*
  * Maximum bounding-box size for naive algorithm
  */
-#define MAX_BOX_SIZE 10000
+#define MAX_BOX_SIZE 1000
 
 /*
  * Maximum number of sets that should be considered

multirow/include/multirow/util.h

@@ -78,6 +78,7 @@ void time_printf(const char *fmt,
 #define UNUSED(x) (void)(x)
 
 #define min(x,y) ((x) < (y) ? (x) : (y))
+#define max(x,y) ((x) > (y) ? (x) : (y))
 #define sign(x) ((x) > 0 ? 1 : ((x) < 0 ? -1 : 0))
 #define swap(x, y, T) do { T SWAP = x; x = y; y = SWAP; } while (0)
 

multirow/src/cg.c

@@ -945,7 +945,10 @@ int CG_add_multirow_cuts(struct CG *cg,
                 LP_free_row(&cut);
                 goto NEXT_COMBINATION;
             }
-            else abort_iff(rval, "generate failed (cut %d)", count);
+            else if(rval) {
+                dump_tableau(&tableau, count);
+                abort_iff(1, "generate failed (cut %d)", count);
+            }
 
             if_verbose_level dump_cut(&cut, count);
 
@@ -1115,4 +1118,32 @@ int CG_total_nz(const struct Tableau *tableau)
     return total_nz;
 }
 
+int CG_print_model(const struct MultiRowModel *model)
+{
+    int rval = 0;
+    int nrows = model->nrows;
+
+    time_printf("       f = [");
+    for (int i = 0; i < nrows; i++)
+        printf("%20.12lf ", model->f[i]);
+    printf("]\n");
+
+    for (int i = 0; i < model->rays.nrays; i++)
+    {
+        double *ray = LFREE_get_ray(&model->rays, i);
+        double norm = 0;
+
+        time_printf("ray[%3d] = [", i);
+        for (int j = 0; j < nrows; j++)
+        {
+            printf("%20.12lf ", ray[j]);
+            norm += fabs(ray[j]);
+        }
+        printf("] norm=%20.12lf \n", norm);
+    }
+
+CLEANUP:
+    return rval;
+}
+
 #endif // TEST_SOURCE

multirow/src/lfree2d.c

@@ -559,6 +559,41 @@ CLEANUP:
     return rval;
 }
 
+int LFREE_print_set(const struct ConvLFreeSet *set)
+{
+    int rval = 0;
+
+    log_debug("f=%12.6lf %12.6lf\n", set->f[0], set->f[1]);
+
+    for (int i = 0; i < set->rays.nrays; i++)
+    {
+        double *ray = LFREE_get_ray(&set->rays, i);
+        log_debug("ray[%-3d] = ", i);
+        for (int j = 0; j < set->nrows; j++)
+            printf("%20.12lf ", ray[j]);
+        printf("\n");
+    }
+
+    for (int i = 0; i < set->rays.nrays; i++)
+    {
+        double beta = set->beta[i];
+        log_debug("beta[%-3d] = %20.12lf\n", i, beta);
+    }
+
+    for (int i = 0; i < set->rays.nrays; i++)
+    {
+        double *ray = LFREE_get_ray(&set->rays, i);
+        double beta = set->beta[i];
+        log_debug("vertex[%-3d] = ", i);
+        for (int j = 0; j < set->nrows; j++)
+            printf("%20.12lf ", ray[j] * beta + set->f[j]);
+        printf("\n");
+    }
+
+CLEANUP:
+    return rval;
+}
+
 void LFREE_free_conv(struct ConvLFreeSet *lfree)
 {
     if(!lfree) return;

multirow/src/linalg.c

@@ -0,0 +1,52 @@
+/*
+ * Copyright (C) 2016 Álinson Santos Xavier <isoron@gmail.com>
+ *
+ * This file is part of Loop Habit Tracker.
+ *
+ * Loop Habit Tracker is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by the
+ * Free Software Foundation, either version 3 of the License, or (at your
+ * option) any later version.
+ *
+ * Loop Habit Tracker is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
+ * more details.
+ *
+ * You should have received a copy of the GNU General Public License along
+ * with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <math.h>
+#include <cblas.h>
+#include <lapacke.h>
+#include <multirow/util.h>
+
+void LINALG_scale(int n, double *x, double alpha)
+{
+    cblas_dscal(n, alpha, x, 1);
+}
+
+double LINALG_dot(int n, double *x, double *y)
+{
+    return cblas_ddot(n, x, 1, y, 1);
+}
+
+double LINALG_norm(int n, double *x)
+{
+    return cblas_dasum(n, x, 1);
+}
+
+int LINALG_solve(int n, int m, double *A, double *b, double *x)
+{
+    int rval = 0;
+    double A_copy[n * m];
+    memcpy(x, b, m * sizeof(double));
+    memcpy(A_copy, A, n * m * sizeof(double));
+
+    LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', m, n, 1, A_copy, n, x, 1);
+    abort_if(rval, "LAPACKE_dgesv failed");
+
+CLEANUP:
+    return rval;
+}

multirow/src/lp.c

@@ -35,6 +35,7 @@ int LP_open(struct LP *lp)
 
     CPXsetintparam(lp->cplex_env, CPX_PARAM_DATACHECK, CPX_ON);
     CPXsetintparam(lp->cplex_env, CPX_PARAM_NUMERICALEMPHASIS, CPX_ON);
+    CPXsetintparam(lp->cplex_env, CPX_PARAM_THREADS, 1);
     CPXsetlogfile(lp->cplex_env, 0);
 
 CLEANUP:

multirow/tests/linalg-test.cpp

@@ -0,0 +1,107 @@
+/* Copyright (c) 2015-2017 Alinson Xavier
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gtest/gtest.h>
+
+extern "C" {
+#include <multirow/util.h>
+#include <multirow/linalg.h>
+}
+
+const double E = 1e-6;
+
+TEST(LinAlgTest, scale_test)
+{
+    double x[] = {1.0, 2.0, 3.0};
+    LINALG_scale(3, x, 2.0);
+    EXPECT_NEAR(x[0], 2.0, E);
+    EXPECT_NEAR(x[1], 4.0, E);
+    EXPECT_NEAR(x[2], 6.0, E);
+}
+
+TEST(LinAlgTest, dot_test)
+{
+    double x[] = { 1.0, 2.0, 3.0 };
+    double y[] = { 3.0, 4.0, 5.0 };
+    double dot = LINALG_dot(3, x, y);
+    EXPECT_NEAR(dot, 26.0, E);
+}
+
+TEST(LinAlgTest, norm_test)
+{
+    double x[] = { 1.0,  2.0, -3.0 };
+    double y[] = { 3.0, -4.0, -5.0 };
+    double x_norm = LINALG_norm(3, x);
+    double y_norm = LINALG_norm(3, y);
+    EXPECT_NEAR(x_norm, 6.0, E);
+    EXPECT_NEAR(y_norm, 12.0, E);
+}
+
+TEST(LinAlgTest, solve_test)
+{
+    int rval = 0;
+
+    double A[] = {
+            2.0, 1.0,  3.0,
+            2.0, 6.0,  8.0,
+            6.0, 8.0, 18.0,
+    };
+    double b[] = { 1.0, 3.0, 5.0 };
+    double x[] = { 0.0, 0.0, 0.0 };
+
+    rval = LINALG_solve(3, 3, A, b, x);
+    abort_if(rval, "LINALG_solve failed");
+
+    // Should compute x correctly
+    EXPECT_NEAR(x[0], 0.3, E);
+    EXPECT_NEAR(x[1], 0.4, E);
+    EXPECT_NEAR(x[2], 0.0, E);
+
+    // Should not modify A and b
+    EXPECT_EQ(A[0], 2.0);
+    EXPECT_EQ(A[1], 1.0);
+    EXPECT_EQ(A[2], 3.0);
+    EXPECT_EQ(b[0], 1.0);
+    EXPECT_EQ(b[1], 3.0);
+    EXPECT_EQ(b[2], 5.0);
+
+CLEANUP:
+    if(rval) FAIL();
+}
+
+TEST(LinAlgTest, solve_test_2)
+{
+    int rval = 0;
+
+    double A[] = {
+            1.0, 3.0,
+            0.0, 0.0,
+            2.0, 2.0,
+            3.0, 5.0,
+    };
+    double b[] = { 5/2.0, 0.0, 8/3.0, 31/6.0 };
+    double x[2];
+
+    rval = LINALG_solve(2, 4, A, b, x);
+    abort_if(rval, "LINALG_solve failed");
+
+    // Should compute x correctly
+    EXPECT_NEAR(x[0], 3/4.0, E);
+    EXPECT_NEAR(x[1], 7/12.0, E);
+
+CLEANUP:
+    if(rval) FAIL();
+}

onerow/benchmark/CMakeLists.txt

@@ -1,2 +1,9 @@
 add_executable(onerow-benchmark.run src/main.cpp)
-target_link_libraries (onerow-benchmark.run LINK_PUBLIC qxx_static onerow_static m pthread ${GMP_LIBRARIES} ${CPLEX_LIBRARIES})
+target_link_libraries (onerow-benchmark.run LINK_PUBLIC
+    qxx_static
+    onerow_static
+    m
+    pthread
+    ${GMP_LIBRARIES}
+    ${CPLEX_LIBRARIES}
+    ${OpenMP_LIBRARIES})

onerow/benchmark/generate_tables

@@ -3,28 +3,21 @@ for i in out/*yaml; do
     IN=${i/.pre.yaml/}
     IN=${IN/out\//}
     grep -q mip_value $i && continue
-    echo $IN
     grep $IN instances/opt.tab | awk '{ print "mip_value:\n  "$2 }' >> $i
 done
 
-echo "                                     TABLE 1                                    "
+printf "%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n" instance cutsmir cutsw origgap mirperf wperf mircontrib wcontrib wimprov wtime > 'tables/gap.csv'
-echo "--------------------------------------------------------------------------------"
-printf "%s,%s,%s,%s,%s,%s,%s,%s,%s,%s\n" instance cutsmir cutsw origgap mirperf wperf mircontrib wcontrib wimprov wtime | tee 'tables/gap.csv'
 for i in out/*.yaml; do
     IN=${i/.yaml/}
     IN=${IN/out\//}
     printf "%s," $IN
     scripts/gap.rb $i
-done | sort -t',' -nrsk 8 | sed -e 's/,$//g' | tee -a 'tables/gap.csv'
+done | sed -e 's/,$//g' >> 'tables/gap.csv'
 
-echo
+printf "%s,%s,%s,%s,%s,%s\n" instance cutsmir cutsw mirt wedget avgm > 'tables/speed.csv'
-echo
-echo "                                     TABLE 2                                    "
-echo "--------------------------------------------------------------------------------"
-printf "%s,%s,%s,%s,%s,%s\n" instance cutsmir cutsw mirt wedget avgm | tee 'tables/speed.csv'
 for i in out/*.yaml; do
     IN=${i/.yaml/}
     IN=${IN/out\//}
     printf "%s," $IN
     scripts/speed.rb $i
-done | sort | sed -e 's/,$//g' | tee -a 'tables/speed.csv'
+done | sort | sed -e 's/,$//g' >> 'tables/speed.csv'

onerow/benchmark/src/main.cpp

@@ -98,6 +98,8 @@ int main(int argc, char **argv)
 
 	Stats::set_input_filename(string(input_filename));
 
+	time_printf("Using OpenMP (%d threads)\n", omp_get_max_threads());
+
 	// reads input file
 	time_printf("Reading input file: %s...\n", input_filename);
 	status = CPXreadcopyprob(env, lp, input_filename, NULL);
@@ -149,21 +151,21 @@ int main(int argc, char **argv)
 	if (enable_gomory_cuts)
 	{
 		time_printf("Generating Gomory cuts...\n");
-		cplexHelper.add_single_row_cuts<GomoryCutGenerator>();
+		cplexHelper.add_single_row_cuts<GomoryCutGenerator>(0);
 		cplexHelper.solve(true);
 	}
 
 	if (enable_mir_cuts)
 	{
 		time_printf("Generating MIR cuts...\n");
-		cplexHelper.add_single_row_cuts<MIRCutGenerator>();
+		cplexHelper.add_single_row_cuts<MIRCutGenerator>(0);
 		cplexHelper.solve(true);
 	}
 
 	if (enable_wedge_cuts)
 	{
 		time_printf("Generating wedge cuts...\n");
-		cplexHelper.add_single_row_cuts<WedgeCutGenerator>();
+		cplexHelper.add_single_row_cuts<WedgeCutGenerator>(MAX_GOOD_ROWS);
 		cplexHelper.solve(true);
 	}
 

onerow/library/include/onerow/cplex_helper.hpp

@@ -69,10 +69,11 @@ public:
 	 * cuts as possible. The cuts are generated by the provided generator class.
 	 *
 	 * @tparam Generator Class used to generate the cuts.
+	 * @param max_rows The maximum number of rows to consider.
 	 * @returns The number of cuts added.
 	 */
 	template<class Generator>
-	int add_single_row_cuts();
+	int add_single_row_cuts(int max_rows);
 
 	/**
 	 * Gets a single row from the current tableau.
@@ -90,7 +91,7 @@ public:
 	 *
 	 * @returns The solution status, as returned by CPXgetstat.
 	 */
-	int solve(bool save_stats = false);
+    void solve(bool save_stats = false);
 
 	void dump_constraint(const Constraint &c, const char *msg = "");
 
@@ -107,7 +108,7 @@ public:
 
 	void print_solution(double *x);
 
-	void find_good_rows();
+	void find_good_rows(int max_rows);
 
 	int n_rows;
 	int n_cols;
@@ -134,6 +135,9 @@ public:
 
 	int n_good_rows;
 	int *good_rows;
+
+	double *reduced_costs;
+	double cost_cutoff;
 };
 
 #include "cplex_helper.tpp"

onerow/library/include/onerow/cplex_helper.tpp

@@ -30,12 +30,17 @@ using std::endl;
 
 
 template<class Generator>
-int CplexHelper::add_single_row_cuts()
+int CplexHelper::add_single_row_cuts(int max_rows)
 {
 	total_cuts = 0;
 
-	if(n_good_rows < 0)
+	if(n_good_rows > 0)
-		find_good_rows();
+	{
+		n_good_rows = 0;
+		delete good_rows;
+	}
+
+	find_good_rows(max_rows);
 
 	eta_reset();
 	eta_count = 0;
@@ -48,7 +53,6 @@ int CplexHelper::add_single_row_cuts()
 	for (int i = 0; i < n_good_rows; i++)
 	{
 		Row *row = get_tableau_row(good_rows[i]);
-		//Row *row = good_rows[i];
 
 		Generator generator(*row);
 

onerow/library/include/onerow/params.hpp

@@ -25,16 +25,14 @@ const long REDUCE_FACTOR_RHS          = 1000000;
 const long REDUCE_FACTOR_R1           = 1000;
 const long REDUCE_FACTOR_COEFFICIENT  = 1000000;
 
-const int MAX_CUT_DEPTH = 200;
+const int MAX_R1_RAYS   = 1000000;
+const int MAX_CUT_DEPTH = 1000000;
+const int MAX_GOOD_ROWS = 1000000;
 
 const int ETA_UPDATE_INTERVAL = 300;
 const unsigned int MAX_CUT_BUFFER_SIZE = 100;
 
-const double COEFFICIENT_SAFETY_MARGIN = 0.00001;
-
 #define INTERSECTION_CUT_USE_DOUBLE
-
-#define ENABLE_TRIVIAL_LIFTING
 // #define ENABLE_EXTENDED_STATISTICS
 // #define PRETEND_TO_ADD_CUTS
 

onerow/library/include/onerow/single_row_cut_generator.hpp

@@ -66,6 +66,8 @@ struct Row {
 	int basic_var_index;
 
 	bool* is_integer;
+	double* reduced_costs;
+	double cost_cutoff;
 };
 
 /**

onerow/library/src/cplex_helper.cpp

@@ -29,6 +29,7 @@
 #include <onerow/geometry.hpp>
 #include <onerow/stats.hpp>
 #include <onerow/params.hpp>
+#include <cstring>
 
 using std::cout;
 using std::endl;
@@ -39,7 +40,7 @@ static bool debug = false;
 CplexHelper::CplexHelper(CPXENVptr _env, CPXLPptr _lp) :
 		env(_env), lp(_lp), is_integer(0), n_cuts(0), n_rows(0), ub(0), lb(0),
 		cstat(0), cplex_rows(0), first_solution(0), current_solution(0),
-		optimal_solution(0), current_round(0), n_good_rows(-1)
+		optimal_solution(0), current_round(0), n_good_rows(-1), reduced_costs(0)
 {
 }
 
@@ -260,6 +261,9 @@ Row* CplexHelper::get_tableau_row(int index)
 	row->is_integer = is_integer;
 	row->c = cplex_row_to_constraint(cplex_rows[index]);
 
+	row->reduced_costs = reduced_costs;
+	row->cost_cutoff = cost_cutoff;
+
 	if (optimal_solution)
 		assert(cplex_rows[index].get_violation(optimal_solution) <= 0.001);
 
@@ -272,7 +276,7 @@ Row* CplexHelper::get_tableau_row(int index)
 	return row;
 }
 
-int CplexHelper::solve(bool should_end_round)
+void CplexHelper::solve(bool should_end_round)
 {
 	// Optimize
 	int status = CPXlpopt(env, lp);
@@ -312,8 +316,6 @@ int CplexHelper::solve(bool should_end_round)
 
 	if (should_end_round)
 		Stats::set_solution(current_round++, objval, string(buffer));
-
-	return objval;
 }
 
 
@@ -364,6 +366,16 @@ void CplexHelper::read_basis()
 	CPXgetlb(env, lp, lb, 0, n_cols - 1);
 	CPXgetbase(env, lp, cstat, rstat);
 
+	reduced_costs = new double[n_cols];
+	CPXgetdj(env, lp, reduced_costs, 0, n_cols-1);
+
+	cost_cutoff = -INFINITY;
+	double *costs_copy = new double[n_cols];
+	memcpy(costs_copy, reduced_costs, sizeof(double) * n_cols);
+	std::sort(costs_copy, costs_copy + n_cols, std::greater<double>());
+	if(n_cols > MAX_R1_RAYS) cost_cutoff = costs_copy[MAX_R1_RAYS];
+	delete costs_copy;
+
 	cplex_rows = new CplexRow[n_rows];
 	assert(cplex_rows != 0);
 
@@ -509,10 +521,10 @@ void CplexHelper::eta_print()
 FINISHED:;
 }
 
-void CplexHelper::find_good_rows()
+void CplexHelper::find_good_rows(int max_rows)
 {
-	n_good_rows = 0;
+	bool *is_good = new bool[n_rows];
-	good_rows = new int[n_rows];
+	double *fractionality = new double[n_rows];
 
 	time_printf("Finding interesting rows...\n");
 
@@ -521,14 +533,37 @@ void CplexHelper::find_good_rows()
 	{
 		Row *row = get_tableau_row(i);
 
-		if (row->c.pi_zero.frac() != 0 &&
+		fractionality[i] = row->c.pi_zero.frac().get_double();
-				row->is_integer[row->basic_var_index])
+		fractionality[i] = fabs(fractionality[i] - 0.5);
-		{
+
-			good_rows[n_good_rows++] = i;
+		is_good[i] = true;
-		}
+		if (row->c.pi_zero.frac() == 0) is_good[i] = false;
+		if (!row->is_integer[row->basic_var_index]) is_good[i] = false;
 
 		delete row;
 	}
 
-	time_printf("	 %d rows found\n", n_good_rows, n_rows);
+	if(max_rows > 0)
+	{
+		double frac_cutoff = 1.0;
+		std::sort(fractionality, fractionality + n_rows);
+		if (n_rows > max_rows) frac_cutoff = fractionality[max_rows];
+
+		for (int i = 0; i < n_rows; i++)
+			if (fractionality[i] > frac_cutoff)
+				is_good[i] = false;
+	}
+
+	n_good_rows = 0;
+	good_rows = new int[n_rows];
+	for (int i = 0; i < n_rows; i++)
+	{
+		if (!is_good[i]) continue;
+		good_rows[n_good_rows++] = i;
+		if(max_rows > 0 && n_good_rows >= max_rows) break;
+	}
+
+	delete is_good;
+	delete fractionality;
+	time_printf("	 %d rows found\n", n_good_rows);
 }

onerow/library/src/wedge_cut_generator.cpp

@@ -296,6 +296,9 @@ void WedgeCutGenerator::eval_next()
 		}
 
 		int r1_index = row.c.pi.index(r1_offset);
+		if(row.reduced_costs[r1_index] < row.cost_cutoff)
+			continue;
+
 		r1[0] = -row.c.pi.value(r1_offset).reduce(REDUCE_FACTOR_R1);
 		r1[1] = 1;