Implement Espinoza's algorithm; other minor improvements

CMakeLists.txt

@@ -9,8 +9,8 @@ find_package(GMP REQUIRED)
 find_package(OpenMP REQUIRED)
 
 set(CMAKE_CXX_STANDARD 11)
-set(CMAKE_CXX_FLAGS  "${CMAKE_CXX_FLAGS} -Werror ${OpenMP_CXX_FLAGS}")
-set(CMAKE_C_FLAGS  "${CMAKE_CXX_FLAGS} -Werror ${OpenMP_C_FLAGS}")
+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)

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,9 @@ 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
+SAMPLES_SLOW=10
+SAMPLES_MEDIUM=100
+SAMPLES_FAST=1000
 SEED=1240
 
 # ORIGINAL
@@ -36,8 +35,11 @@ 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 $SAMPLES_FAST --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 $SAMPLES_FAST --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
@@ -47,11 +49,17 @@ $RUN $COMMON_OPTS --samples $SAMPLES_MEDIUM --naive --log $DIR/naive-bbox.log --
 
 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 $SAMPLES_SLOW --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 $SAMPLES_SLOW --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
 
+title MIP Pre-processing
+$RUN $COMMON_OPTS --samples $SAMPLES_SLOW --mip --preprocess --log $DIR/mip-pre.log --stats $DIR/mip-pre.yaml || exit
+
 # SHEAR
 # ------------------------------------------------------------------------------
 
@@ -69,7 +77,7 @@ title Bound Pre-processing + Shear
 $RUN $COMMON_OPTS --samples $SAMPLES_FAST --bound --preprocess --log $DIR/bound-pre.log --stats $DIR/bound-pre.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 $SAMPLES_MEDIUM --bound --log $DIR/bound.log --stats $DIR/bound.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

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,
-                        &value);
+                rval = LIFTING_2D_bound(set->n_halfspaces, set->halfspaces,
+                        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;
-    if(SELECT_MIP_ALGORITHM) algorithm = ALGORITHM_MIP;
-    else if(SELECT_NAIVE_ALGORITHM) algorithm = ALGORITHM_NAIVE;
-
-    if(algorithm == ALGORITHM_NAIVE)
+    int algorithm = -1;
+    
+    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,
-Average (ms),0.565,0.059,141.938,9.341,23.278,
-Median (ms),0.068,0.060,129.200,2.120,23.080,
-Maximum (ms),33.264,0.100,947.600,1769.880,30.080,
-Failure Rate,0.0 \%,0.0 \%,0.0 \%,0.0 \%,0.3 \%,
-Best,37.2 \%,83.6 \%,0.0 \%,0.0 \%,0.0 \%,
-Avg Ratio to Best,9.483,1.015,2449.938,159.449,402.950,
+metric,boundpre,bound,heur,mippre,mip,naivebboxpre,naivebbox,naivefixed50,
+Average (ms),0.078,0.937,0.004,70.080,81.738,6.353,12.490,34.874,
+Median (ms),0.079,0.099,0.004,67.590,73.889,2.900,3.120,34.195,
+Maximum (ms),0.152,55.998,0.008,175.773,1705.441,96.585,2483.652,122.881,
+Failure Rate,0.0 \%,0.0 \%,100.0 \%,0.0 \%,0.0 \%,0.0 \%,0.1 \%,0.3 \%,
+Best,0.0 \%,0.0 \%,100.0 \%,0.0 \%,0.0 \%,0.0 \%,0.0 \%,0.0 \%,
+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,
-Average (ms),4.709,0.062,1906.081,5473.688,23.374,
-Median (ms),0.440,0.064,582.800,21.200,23.160,
-Maximum (ms),485.760,0.108,94620.000,392259.600,32.120,
-Failure Rate,0.0 \%,0.0 \%,1.7 \%,0.1 \%,10.0 \%,
-Best,0.0 \%,100.0 \%,0.0 \%,0.0 \%,0.0 \%,
-Avg Ratio to Best,74.602,1.000,31250.644,87472.010,380.532,
+metric,boundpre,bound,mip,naivebbox,naivefixed50,
+Average (ms),0.080,7.723,659.486,478.703,34.680,
+Median (ms),0.080,0.730,131.680,31.395,34.105,
+Maximum (ms),0.352,800.768,53099.128,9450.563,48.423,
+Failure Rate,0.0 \%,0.0 \%,2.1 \%,0.5 \%,10.0 \%,
+Best,100.0 \%,0.0 \%,0.0 \%,0.0 \%,0.0 \%,
+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/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/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:
@@ -736,4 +737,4 @@ int LP_init_row(struct Row *row, int nz_capacity)
 
     CLEANUP:
     return rval;
-}
+}