diff --git a/infinity/benchmark/2row-cont/gt2.pre.log b/infinity/benchmark/2row-cont/gt2.pre.log index b3da870..92ba46c 100644 --- a/infinity/benchmark/2row-cont/gt2.pre.log +++ b/infinity/benchmark/2row-cont/gt2.pre.log @@ -1,30 +1,25 @@ -[ 0.00] multirow (git-4d20687d04e75dc63c7bb84cbbfeb2261c6a6ab2) -[ 0.00] mf-hpc-1-7.local 2017-03-11 22:40 -[ 0.00] Compile-time parameters: -[ 0.00] EPSILON: 1.000000e-08 -[ 0.00] GREEDY_BIG_E: 1.000000e+03 -[ 0.00] GREEDY_MAX_GAP: 1.000000e-04 -[ 0.00] Command line arguments: -[ 0.00] ../build/infinity.run --mir --greedy --problem instances/gt2.pre.mps.gz --basis bases/gt2.pre.bas --log 2row-cont/gt2.pre.log --stats 2row-cont/gt2.pre.yaml --solution solutions/gt2.pre.x -[ 0.00] 28 rows, 201 cols -[ 0.00] Storing column types... -[ 0.00] Relaxing integrality... -[ 0.00] Disabling presolve... -[ 0.00] Optimizing... -[ 0.00] opt = 20146.761297 -[ 0.00] Reading tableau rows... -[ 0.00] Adding MIR cuts... -[ 0.00] opt = 20545.739873 -[ 0.00] opt = 20565.114738 -[ 0.00] opt = 20628.064931 -[ 0.00] opt = 20756.821715 -[ 0.00] Optimizing... -[ 0.00] opt = 20756.821715 -[ 0.00] Adding greedy intersection cuts (2 rows)... -[ 0.00] Finding combinations... -[ 0.00] 249 combinations [0.05] -[ 0.00] opt = 20815.431583 -[ 0.05] opt = 21061.100847 -[ 0.26] Optimizing... -[ 0.26] opt = 21061.100847 -[ 0.26] Writing stats to file 2row-cont/gt2.pre.yaml... +[ 0.04] Reading problem instances/gt2.pre.mps.gz... +[ 0.41] 28 rows, 201 cols +[ 0.41] Storing column types... +[ 0.42] Relaxing integrality... +[ 0.43] Disabling presolve... +[ 0.43] Reading basis from file bases/gt2.pre.bas... +[ 0.47] Optimizing... +[ 0.74] opt = 20146.761297 +[ 0.74] Reading tableau rows... +[ 0.80] Reading solution from solutions/gt2.pre.x +[ 0.83] Adding MIR cuts... +[ 1.09] opt = 20545.739873 +[ 1.17] opt = 20565.114738 +[ 1.21] opt = 20628.064931 +[ 1.28] opt = 20756.821715 +[ 1.28] Optimizing... +[ 1.29] opt = 20756.821715 +[ 1.29] Adding greedy intersection cuts (2 rows)... +[ 1.29] Finding combinations... +[ 1.31] 249 combinations [0.05] +[ 1.75] opt = 20815.431583 +[ 3.27] opt = 21061.100847 +[ 7.74] Optimizing... +[ 7.74] opt = 21061.100847 +[ 7.74] Writing stats to file 2row-cont/gt2.pre.yaml... diff --git a/infinity/benchmark/2row-cont/gt2.pre.yaml b/infinity/benchmark/2row-cont/gt2.pre.yaml index 821e5cb..bdb8c02 100644 --- a/infinity/benchmark/2row-cont/gt2.pre.yaml +++ b/infinity/benchmark/2row-cont/gt2.pre.yaml @@ -13,9 +13,9 @@ added-cuts: 1: 7 2: 2 user-cpu-time: - 0: 0.000 - 1: 0.000 - 2: 0.267 + 0: 0.871 + 1: 0.548 + 2: 6.451 time-per-cut: - 1: 0.000 - 2: 0.001 + 1: 0.037 + 2: 0.026 diff --git a/infinity/benchmark/2row-cont/harp2.pre.log b/infinity/benchmark/2row-cont/harp2.pre.log index f43ca6b..0309caa 100644 --- a/infinity/benchmark/2row-cont/harp2.pre.log +++ b/infinity/benchmark/2row-cont/harp2.pre.log @@ -1,40 +1,17 @@ -[ 0.00] multirow (git-4d20687d04e75dc63c7bb84cbbfeb2261c6a6ab2) -[ 0.00] mf-hpc-1-7.local 2017-03-11 22:40 -[ 0.00] Compile-time parameters: -[ 0.00] EPSILON: 1.000000e-08 -[ 0.00] GREEDY_BIG_E: 1.000000e+03 -[ 0.00] GREEDY_MAX_GAP: 1.000000e-04 -[ 0.00] Command line arguments: -[ 0.00] ../build/infinity.run --mir --greedy --problem instances/harp2.pre.mps.gz --basis bases/harp2.pre.bas --log 2row-cont/harp2.pre.log --stats 2row-cont/harp2.pre.yaml --solution solutions/harp2.pre.x +[ 0.00] Reading problem instances/harp2.pre.mps.gz... [ 0.00] 92 rows, 1049 cols [ 0.00] Storing column types... [ 0.00] Relaxing integrality... [ 0.00] Disabling presolve... -[ 0.01] Optimizing... -[ 0.01] opt = -74325169.345138 -[ 0.01] Reading tableau rows... +[ 0.00] Reading basis from file bases/harp2.pre.bas... +[ 0.00] Optimizing... +[ 0.00] opt = -74325169.345138 +[ 0.00] Reading tableau rows... +[ 0.01] Reading solution from solutions/harp2.pre.x [ 0.01] Adding MIR cuts... [ 0.01] opt = -74315010.901620 -[ 0.01] opt = -74314964.934360 -[ 0.01] opt = -74313681.237627 -[ 0.01] opt = -74313582.579756 [ 0.01] opt = -74293298.124331 -[ 0.01] opt = -74292514.336406 -[ 0.01] opt = -74281276.834560 -[ 0.01] opt = -74278931.707997 -[ 0.01] opt = -74274793.244156 -[ 0.01] opt = -74268371.376509 -[ 0.01] opt = -74264893.739723 -[ 0.01] opt = -74240861.206190 -[ 0.01] opt = -74222594.677631 -[ 0.01] Optimizing... -[ 0.01] opt = -74222594.677631 -[ 0.01] Adding greedy intersection cuts (2 rows)... -[ 0.01] Finding combinations... -[ 0.02] 483 combinations [0.05] -[ 0.15] opt = -74222542.959159 -[ 0.45] opt = -74222502.242291 -[ 1.44] opt = -74222380.900680 -[ 2.56] Optimizing... -[ 2.56] opt = -74222380.900680 -[ 2.56] Writing stats to file 2row-cont/harp2.pre.yaml... +[ 0.01] Cut cuts off known integral solution: -0.26968016 >= -0.42512401 (0) (/Users/axavier/.shared/research/multirow/multirow/src/cg.c:319) +[ 0.01] check_cut failed (1) (/Users/axavier/.shared/research/multirow/multirow/src/cg.c:339) +[ 0.01] add_cut failed (1) (/Users/axavier/.shared/research/multirow/multirow/src/cg.c:644) +[ 0.01] CG_add_single_row_cuts failed (1) (/Users/axavier/.shared/research/multirow/infinity/benchmark/src/main.c:397) diff --git a/infinity/benchmark/src/main.c b/infinity/benchmark/src/main.c index f76eb22..b69c88a 100644 --- a/infinity/benchmark/src/main.c +++ b/infinity/benchmark/src/main.c @@ -25,14 +25,14 @@ #include #include -#include +#include int ENABLE_LIFTING = 0; int MIN_N_ROWS = 2; int MAX_N_ROWS = 2; -int DUMP_CUT = 0; +int SHOULD_DUMP_CUTS = 0; int DUMP_CUT_N = 0; int GENERATE_MIR = 0; @@ -415,8 +415,8 @@ int main(int argc, { log_info("Adding greedy intersection cuts (%d rows)...\n", k); - rval = CG_add_multirow_cuts(cg, k, (MultirowGeneratorCallback) - GREEDY_generate_cut); + rval = CG_add_multirow_cuts(cg, k, + (MultiRowGeneratorCallback) INFINITY_generate_cut); abort_if(rval, "CG_add_multirow_cuts failed"); } diff --git a/infinity/library/CMakeLists.txt b/infinity/library/CMakeLists.txt index 82e51b5..38cb137 100644 --- a/infinity/library/CMakeLists.txt +++ b/infinity/library/CMakeLists.txt @@ -1,20 +1,22 @@ set(COMMON_SOURCES - src/greedy-2d.c + src/infinity-2d.c src/greedy-nd.c src/greedy-bsearch.c - src/greedy.c - include/infinity/greedy-2d.h + src/infinity.c + include/infinity/infinity-2d.h include/infinity/greedy-nd.h include/infinity/greedy-bsearch.h - include/infinity/greedy.h) + include/infinity/infinity.h) set(TEST_SOURCES - tests/greedy-2d-test.cpp - tests/greedy-nd-test.cpp) + tests/infinity-2d-test.cpp + tests/greedy-nd-test.cpp + tests/infinity-test.cpp) add_library(infinity_static ${COMMON_SOURCES}) set_target_properties(infinity_static PROPERTIES OUTPUT_NAME infinity) target_include_directories (infinity_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) add_executable(infinity-test.run ${COMMON_SOURCES} ${TEST_SOURCES}) +target_compile_options(infinity-test.run PRIVATE "-fpermissive") target_link_libraries(infinity-test.run gtest_main multirow_static lifting_static) diff --git a/infinity/library/include/infinity/greedy.h b/infinity/library/include/infinity/infinity-2d.h similarity index 57% rename from infinity/library/include/infinity/greedy.h rename to infinity/library/include/infinity/infinity-2d.h index 9e02cf1..ee4bc96 100644 --- a/infinity/library/include/infinity/greedy.h +++ b/infinity/library/include/infinity/infinity-2d.h @@ -13,19 +13,12 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef MULTIROW_GREEDY_H -#define MULTIROW_GREEDY_H -int GREEDY_write_sage_file(int nrows, - int nrays, - const double *f, - const double *rays, - const double *bounds, - const char *filename); +#ifndef MULTIROW_INFINITY_2D_H +#define MULTIROW_INFINITY_2D_H -int GREEDY_generate_cut(int nrows, - struct Row **rows, - const char *column_types, - struct Row *cut); +#include -#endif //MULTIROW_GREEDY_H +int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds); + +#endif //MULTIROW_INFINITY_2D_H diff --git a/infinity/library/include/infinity/greedy-2d.h b/infinity/library/include/infinity/infinity.h similarity index 52% rename from infinity/library/include/infinity/greedy-2d.h rename to infinity/library/include/infinity/infinity.h index 267fd6e..1b08782 100644 --- a/infinity/library/include/infinity/greedy-2d.h +++ b/infinity/library/include/infinity/infinity.h @@ -13,25 +13,12 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ +#ifndef MULTIROW_INFINITY_H +#define MULTIROW_INFINITY_H -#ifndef MULTIROW_GREEDY_2D_H -#define MULTIROW_GREEDY_2D_H +#include +#include -int GREEDY_2D_bound(const double *rays, - const double *bounds, - int nrays, - const double *f, - const double *p, - double *epsilon, - double *v1, - double *v2, - int *index1, - int *index2); +int INFINITY_generate_cut(struct Tableau *tableau, struct Row *cut); -int GREEDY_2D_generate_cut(const double *rays, - int nrays, - const double *f, - double *bounds); - - -#endif //MULTIROW_GREEDY_2D_H +#endif //MULTIROW_INFINITY_H diff --git a/infinity/library/src/greedy.c b/infinity/library/src/greedy.c deleted file mode 100644 index abe60b4..0000000 --- a/infinity/library/src/greedy.c +++ /dev/null @@ -1,477 +0,0 @@ -/* Copyright (c) 2015 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 . - */ - -#include -#include - -#include -#include -#include - -#include -#include -#include - -struct SortPair -{ - int index; - void *data; -}; - -static int _qsort_cmp_rays_angle(const void *p1, - const void *p2) -{ - double *r1 = (double *) (((struct SortPair *) p1)->data); - double *r2 = (double *) (((struct SortPair *) p2)->data); - return sign(atan2(r1[0], r1[1]) - atan2(r2[0], r2[1])); -} - -static int sort_rays_angle(double *rays, - int nrays, - double *beta) -{ - int rval = 0; - - int *map = 0; - double *rays_copy = 0; - double *beta_copy = 0; - struct SortPair *pairs = 0; - - pairs = (struct SortPair *) malloc(nrays * sizeof(struct SortPair)); - rays_copy = (double *) malloc(2 * nrays * sizeof(double)); - beta_copy = (double*) malloc(nrays * sizeof(double)); - abort_if(!pairs, "could not allocate pairs"); - abort_if(!rays_copy, "could not allocate rays_copy"); - abort_if(!beta_copy, "could not allocate beta_copy"); - - memcpy(rays_copy, rays, 2 * nrays * sizeof(double)); - memcpy(beta_copy, beta, nrays * sizeof(double)); - - for (int i = 0; i < nrays; i++) - { - pairs[i].index = i; - pairs[i].data = &rays[2 * i]; - } - - qsort(pairs, nrays, sizeof(struct SortPair), _qsort_cmp_rays_angle); - - for (int i = 0; i < nrays; i++) - { - beta[i] = beta_copy[pairs[i].index]; - memcpy(&rays[2 * i], &rays_copy[2 * pairs[i].index], 2 * - sizeof(double)); - } - -CLEANUP: - if (pairs) free(pairs); - if (rays_copy) free(rays_copy); - if (beta_copy) free(beta_copy); - return rval; -} - - -static int scale_rays(double *rays, - int nrays, - double *scale) -{ - int rval = 0; - - for (int i = 0; i < nrays; i++) - { - rays[2*i] *= scale[i]; - rays[2*i+1] *= scale[i]; - } - -CLEANUP: - return rval; -} - - -static void print_row(const struct Row *row) -{ - time_printf("Row:\n"); - for (int i = 0; i < row->nz; i++) - time_printf(" %.4lfx%d\n", row->pi[i], row->indices[i]); - time_printf(" <= %.4lf [%d]\n", row->pi_zero, row->head); -} - -static void print_rays(const int *map, - const int *indices, - const double *f, - const double *rays, - const double *scale, - int nrays, - int nz, - int nrows) -{ - time_printf("Mapping:\n"); - for (int i = 0; i < nz; i++) - time_printf(" %4d: %4d x%-4d (scale=%.4lf)\n", i, map[i], indices[i], - scale[i]); - - time_printf("Origin:\n"); - for (int i = 0; i < nrows; i++) - time_printf(" %20.12lf\n", f[i]); - - time_printf("Rays:\n"); - for (int i = 0; i < nrays; i++) - { - time_printf(" "); - - for (int j = 0; j < nrows; j++) - printf("%20.12lf ", rays[i * nrows + j]); - - printf("angle=%.4lf ", atan2(rays[i * nrows], rays[i * nrows + 1])); - printf("norm=%.4lf ", fabs(rays[i * nrows]) + fabs(rays[i * nrows + 1])); - - printf("[ "); - for (int j = 0; j < nz; j++) - if (map[j] == i) printf("%d ", indices[j]); - printf("]\n"); - } -} - -static void print_cut(const struct Row *cut) -{ - time_printf("Generated cut:\n"); - for (int i = 0; i < cut->nz; i++) - time_printf(" %.4lfx%d\n", cut->pi[i], cut->indices[i]); - time_printf(" <= %.4lf\n", cut->pi_zero); -} - -int GREEDY_write_sage_file(int nrows, - int nrays, - const double *f, - const double *rays, - const double *beta, - const char *filename) -{ - int rval = 0; - - FILE *fsage = fopen(filename, "w"); - abort_iff(!fsage, "could not open %s", filename); - - fprintf(fsage, "f=vector(["); - for (int i = 0; i < nrows; i++) - fprintf(fsage, "%.20lf,", f[i]); - fprintf(fsage, "])\n"); - - fprintf(fsage, "R=matrix([\n"); - for (int i = 0; i < nrays; i++) - { - fprintf(fsage, " ["); - for (int j = 0; j < nrows; j++) - { - fprintf(fsage, "%.20lf,", rays[i * nrows + j]); - } - fprintf(fsage, "],\n"); - } - fprintf(fsage, "])\n"); - - fprintf(fsage, "pi=vector([\n"); - for (int k = 0; k < nrays; k++) - fprintf(fsage, " %.12lf,\n", 1 / beta[k]); - fprintf(fsage, "])\n"); - -CLEANUP: - fclose(fsage); - return rval; -} - -static int create_cut(const int nrows, - const char *column_types, - const int nrays, - const double *f, - const int lfree_nrays, - const double *lfree_rays, - const double *beta, - const double *extracted_rays, - const int nvars, - const int *variable_to_ray, - const int *indices, - const double *scale, - struct Row *cut) -{ - int rval = 0; - - cut->nz = nvars; - cut->pi = (double *) malloc(nvars * sizeof(double)); - cut->indices = (int *) malloc(nvars * sizeof(int)); - abort_if(!cut->pi, "could not allocate cut->pi"); - abort_if(!cut->indices, "could not allocate cut->indices"); - - struct LP lp; - - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); - - rval = GREEDY_create_psi_lp(nrows, lfree_nrays, f, lfree_rays, beta, &lp); - abort_if(rval, "create_psi_lp failed"); - - for (int i = 0; i < nvars; i++) - { - double value; - const double *q = &extracted_rays[variable_to_ray[i] * nrows]; - - if(ENABLE_LIFTING && column_types[indices[i]] == MILP_INTEGER) - { - rval = GREEDY_ND_pi(nrows, lfree_nrays, f, lfree_rays, beta, q, - scale[i], &lp, &value); - abort_if(rval, "GREEDY_ND_pi failed"); - } - else - { - rval = GREEDY_ND_psi(nrows, lfree_nrays, f, lfree_rays, beta, q, - scale[i], &lp, &value); - abort_if(rval, "GREEDY_ND_psi failed"); - } - - log_verbose(" psi[%4d] = %20.12lf %d\n", indices[i], value); - - value *= 1.0001; - value = DOUBLE_max(value, 0.0001); - - cut->indices[i] = indices[i]; - cut->pi[i] = - value; - } - - cut->pi_zero = -1.0; - -CLEANUP: - LP_free(&lp); - return rval; -} - - -#ifndef TEST_SOURCE -int GREEDY_generate_cut(int nrows, - struct Row **rows, - const char *column_types, - struct Row *cut) -{ - int rval = 0; - double *f = 0; - long max_nrays = 0; - - f = (double *) malloc(nrows * sizeof(double)); - abort_if(!f, "could not allocate f"); - - for (int i = 0; i < nrows; i++) - { - f[i] = frac(rows[i]->pi_zero); - if (DOUBLE_eq(f[i], 1.0)) f[i] = 0.0; - - max_nrays += rows[i]->nz; - } - - int nz; - int nrays; - int *variable_to_ray = 0; - int *indices = 0; - double *extracted_rays = 0; - double *scale = 0; - - double *beta = 0; - double *scaled_rays = 0; - - int lfree_nrays = 0; - double *lfree_rays = 0; - - variable_to_ray = (int *) malloc(max_nrays * sizeof(int)); - indices = (int *) malloc(max_nrays * sizeof(int)); - scale = (double *) malloc(max_nrays * sizeof(double)); - extracted_rays = (double *) malloc(nrows * max_nrays * sizeof(double)); - abort_if(!variable_to_ray, "could not allocate variable_to_ray"); - abort_if(!indices, "could not allocate indices"); - abort_if(!scale, "could not allocate scale"); - abort_if(!extracted_rays, "could not allocate extracted_rays"); - - lfree_rays = (double*) malloc(nrows * (max_nrays + 100) * sizeof(double)); - abort_if(!lfree_rays, "could not allocate lfree_rays"); - - rval = CG_extract_rays_from_rows(nrows, rows, extracted_rays, &nrays, - variable_to_ray, scale, indices, &nz); - abort_if(rval, "CG_extract_rays_from_rows failed"); - - for (double norm_cutoff = 0.00; norm_cutoff <= 5.0; norm_cutoff+= 0.1) - { - lfree_nrays = 0; - - for (int i = 0; i < nrays; i++) - { - int keep = 1; - double *r = &extracted_rays[nrows * i]; - - for (int j = 0; j < lfree_nrays; j++) - { - double *q = &extracted_rays[nrows * j]; - double norm = 0; - - for(int k = 0; k < nrows; k++) - norm += fabs(r[k] - q[k]); - - if(norm <= norm_cutoff) - { - keep = 0; - break; - } - } - - if(keep) - { - memcpy(&lfree_rays[nrows * lfree_nrays], r, nrows * sizeof(double)); - lfree_nrays++; - } - } - - log_debug(" norm_cutoff=%8.2lf nrays=%8d\n", norm_cutoff, lfree_nrays); - - if(lfree_nrays < MAX_N_RAYS) break; - } - - if(ENABLE_LIFTING) - { - abort_if(nrows > 3, "not implemented"); - - int n_extra_rays; - double extra_rays[100]; - - if(nrows == 2) - { - extra_rays[0] = 0.0001; - extra_rays[1] = 0.0000; - - extra_rays[2] = -0.0001; - extra_rays[3] = 0.0001; - - extra_rays[4] = -0.0001; - extra_rays[5] = -0.0001; - - n_extra_rays = 3; - } - else if(nrows == 3) - { - extra_rays[0] = 0.0000; - extra_rays[1] = 0.0000; - extra_rays[2] = 0.0001; - - extra_rays[3] = 0.0001; - extra_rays[4] = 0.0000; - extra_rays[5] = -0.0001; - - extra_rays[6] = -0.0001; - extra_rays[7] = 0.0000; - extra_rays[8] = -0.0001; - - extra_rays[ 9] = -0.0001; - extra_rays[10] = -0.0001; - extra_rays[11] = -0.0001; - - n_extra_rays = 4; - } - - - for(int i = 0; i < n_extra_rays; i++) - { - double *r = &extra_rays[nrows * i]; - - double scale; - int found, index; - - rval = CG_find_ray(nrows, lfree_rays, lfree_nrays, r, &found, - &scale, &index); - abort_if(rval, "CG_find_ray failed"); - - if(!found) { - memcpy(&lfree_rays[nrows * lfree_nrays], r, nrows * - sizeof(double)); - lfree_nrays++; - } - } - } - - - if(lfree_nrays < 3) - { - rval = ERR_NO_CUT; - cut->pi = 0; - cut->indices = 0; - goto CLEANUP; - } - - log_debug("Extracted %d rays\n", lfree_nrays); - if_verbose_level - { - print_rays(variable_to_ray, indices, f, extracted_rays, scale, nrays, - nz, nrows); - } - - beta = (double *) malloc((lfree_nrays) * sizeof(double)); - abort_if(!beta, "could not allocate beta"); - - log_verbose("Computing lattice-free set...\n"); - - if(nrows == 2) - { - rval = sort_rays_angle(lfree_rays, lfree_nrays, beta); - abort_if(rval, "sort_rays_angle failed"); - - rval = GREEDY_2D_generate_cut(lfree_rays, lfree_nrays, f, beta); - if(rval) - { - rval = ERR_NO_CUT; - goto CLEANUP; - } - } - else - { - rval = GREEDY_ND_generate_cut(nrows, lfree_nrays, f, lfree_rays, beta); - if(rval) - { - rval = ERR_NO_CUT; - goto CLEANUP; - } - } - - - if(DUMP_CUT) - { - char filename[100]; - sprintf(filename, "cut-%03d.sage", DUMP_CUT_N++); - - time_printf("Writing %s...\n", filename); - rval = GREEDY_write_sage_file(nrows, lfree_nrays, f, lfree_rays, beta, - filename); - abort_if(rval, "GREEDY_write_sage_file failed"); - } - - rval = create_cut(nrows, column_types, nrays, f, lfree_nrays, lfree_rays, - beta, extracted_rays, nz, variable_to_ray, indices, scale, cut); - abort_if(rval, "create_cut failed"); - - if_verbose_level print_cut(cut); - -CLEANUP: - if (f) free(f); - if (variable_to_ray) free(variable_to_ray); - if (beta) free(beta); - if (indices) free(indices); - if (scale) free(scale); - if (scaled_rays) free(scaled_rays); - if (lfree_rays) free(lfree_rays); - return rval; -} -#endif // TEST_SOURCE diff --git a/infinity/library/src/greedy-2d.c b/infinity/library/src/infinity-2d.c similarity index 95% rename from infinity/library/src/greedy-2d.c rename to infinity/library/src/infinity-2d.c index 8a8fd12..af3d2f8 100644 --- a/infinity/library/src/greedy-2d.c +++ b/infinity/library/src/infinity-2d.c @@ -21,9 +21,9 @@ #include #include #include -#include +#include -#include +#include static int get_bounding_box(int nrows, int nrays, @@ -390,18 +390,16 @@ CLEANUP: return rval; } -#ifndef TEST_SOURCE - -int GREEDY_2D_bound(const double *rays, - const double *bounds, - int nrays, - const double *f, - const double *p, - double *epsilon, - double *v1, - double *v2, - int *index1, - int *index2) +static int bound(const double *rays, + const double *bounds, + int nrays, + const double *f, + const double *p, + double *epsilon, + double *v1, + double *v2, + int *index1, + int *index2) { int rval = 0; @@ -518,19 +516,21 @@ int GREEDY_2D_bound(const double *rays, } } - CLEANUP: +CLEANUP: return rval; } -int GREEDY_2D_generate_cut(const double *original_rays, - const int nrays, - const double *f, - double *bounds) +#ifndef TEST_SOURCE + + +int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) { - log_verbose("GREEDY_2D_generate_cut\n"); + log_verbose("INFINITY_2D_generate_cut\n"); int rval = 0; int count = 0; + int nrays = model->nrays; + double *f = model->f; double *scale = 0; double *rays = 0; @@ -545,7 +545,7 @@ int GREEDY_2D_generate_cut(const double *original_rays, abort_if(!rays, "could not allocate rays"); abort_if(!scale, "could not allocate scale"); - memcpy(rays, original_rays, 2 * nrays * sizeof(double)); + memcpy(rays, model->rays, 2 * nrays * sizeof(double)); rval = scale_to_chull(rays, nrays, scale); abort_if(rval, "scale_to_chull failed"); @@ -581,9 +581,8 @@ int GREEDY_2D_generate_cut(const double *original_rays, log_verbose(" p=%.2lf %.2lf\n", p[0], p[1]); - rval = GREEDY_2D_bound(rays, bounds, nrays, f, p, &epsilon, v1, v2, - &i1, &i2); - abort_if(rval, "GREEDY_2D_bound failed"); + rval = bound(rays, bounds, nrays, f, p, &epsilon, v1, v2, &i1, &i2); + abort_if(rval, "bound failed"); log_verbose(" epsilon=%.2lf\n", epsilon); @@ -720,8 +719,6 @@ int GREEDY_2D_generate_cut(const double *original_rays, } if(is_split) break; - - } for(int i=0; i. + */ + +#include +#include + +#include +#include +#include + +#include +#include +#include + +struct SortPair +{ + int index; + void *data; +}; + +/** + * Compares two rays according to their angles. + * + * The input are two SortPairs, where the data points to a double[2]. If the + * two rays are pointing to exactly opposite direction, returns 1. + * + * @param p1 a pointer to a SortPair containing the first ray. + * @param p2 a pointer to a SortPair containing the second ray. + * @return -1, 1 or 0, if the second ray is in clockwise, counter-clockwise + * or aligned, respectively, to the first ray. + */ +static int _qsort_cmp_rays_angle(const void *p1, const void *p2) +{ + double *r1 = (double *) (((struct SortPair *) p1)->data); + double *r2 = (double *) (((struct SortPair *) p2)->data); + return sign(atan2(r1[0], r1[1]) - atan2(r2[0], r2[1])); +} + +/** + * Sorts a list of rays according to their angle. + * + * @param rays the rays to be sorted + * @param nrays the number of rays in the list + * @param beta a list of doubles, to be sorted together with the rays + * + * @return zero if successful, non-zero otherwise + */ +static int sort_rays_angle(double *rays, int nrays, double *beta) +{ + int rval = 0; + + double *rays_copy = 0; + double *beta_copy = 0; + struct SortPair *pairs = 0; + + pairs = (struct SortPair *) malloc(nrays * sizeof(struct SortPair)); + rays_copy = (double *) malloc(2 * nrays * sizeof(double)); + beta_copy = (double *) malloc(nrays * sizeof(double)); + abort_if(!pairs, "could not allocate pairs"); + abort_if(!rays_copy, "could not allocate rays_copy"); + abort_if(!beta_copy, "could not allocate beta_copy"); + + memcpy(rays_copy, rays, 2 * nrays * sizeof(double)); + memcpy(beta_copy, beta, nrays * sizeof(double)); + + for (int i = 0; i < nrays; i++) + { + pairs[i].index = i; + pairs[i].data = &rays[2 * i]; + } + + qsort(pairs, (size_t) nrays, sizeof(struct SortPair), + _qsort_cmp_rays_angle); + + for (int i = 0; i < nrays; i++) + { + beta[i] = beta_copy[pairs[i].index]; + memcpy(&rays[2 * i], &rays_copy[2 * pairs[i].index], + 2 * sizeof(double)); + } + +CLEANUP: + if (pairs) free(pairs); + if (rays_copy) free(rays_copy); + if (beta_copy) free(beta_copy); + return rval; +} + +static void print_row(const struct Row *row) +{ + time_printf("Row:\n"); + for (int i = 0; i < row->nz; i++) + time_printf(" %.4lfx%d\n", row->pi[i], row->indices[i]); + time_printf(" <= %.4lf [%d]\n", row->pi_zero, row->head); +} + +static void print_rays(const struct Tableau *tableau, + const struct RayMap *map, + const double *f, + const double *rays, + int nrays) +{ + int nrows = tableau->nrows; + + time_printf("Ray map:\n"); + for (int i = 0; i < map->nvars; i++) + time_printf(" %4d: %4d x%-4d (scale=%.4lf)\n", i, + map->variable_to_ray[i], map->indices[i], map->ray_scale[i]); + + time_printf("Origin:\n"); + for (int i = 0; i < nrows; i++) + time_printf(" %20.12lf\n", f[i]); + + time_printf("Rays:\n"); + for (int i = 0; i < nrays; i++) + { + time_printf(" "); + + for (int j = 0; j < nrows; j++) + printf("%20.12lf ", rays[i * nrows + j]); + + printf("angle=%.4lf ", atan2(rays[i * nrows], rays[i * nrows + 1])); + printf("norm=%.4lf ", + fabs(rays[i * nrows]) + fabs(rays[i * nrows + 1])); + + printf("[ "); + for (int j = 0; j < map->nvars; j++) + if (map->variable_to_ray[j] == i) + printf("%d ", map->indices[j]); + printf("]\n"); + } +} + +static void print_cut(const struct Row *cut) +{ + time_printf("Generated cut:\n"); + for (int i = 0; i < cut->nz; i++) + time_printf(" %.4lfx%d\n", cut->pi[i], cut->indices[i]); + time_printf(" <= %.4lf\n", cut->pi_zero); +} + +static int create_cut(const struct Tableau *tableau, + const struct RayMap *map, + const double *f, + const int lfree_nrays, + const double *lfree_rays, + const double *beta, + struct Row *cut) +{ + int rval = 0; + int nvars = map->nvars; + int nrows = tableau->nrows; + + cut->nz = nvars; + cut->pi = (double *) malloc(nvars * sizeof(double)); + cut->indices = (int *) malloc(nvars * sizeof(int)); + abort_if(!cut->pi, "could not allocate cut->pi"); + abort_if(!cut->indices, "could not allocate cut->indices"); + + struct LP lp; + + rval = LP_open(&lp); + abort_if(rval, "LP_open failed"); + + rval = GREEDY_create_psi_lp(nrows, lfree_nrays, f, lfree_rays, beta, &lp); + abort_if(rval, "create_psi_lp failed"); + + for (int i = 0; i < nvars; i++) + { + double value; + const double *q = &map->rays[map->variable_to_ray[i] * nrows]; + + if (ENABLE_LIFTING && + tableau->column_types[map->indices[i]] == MILP_INTEGER) + { + rval = GREEDY_ND_pi(nrows, lfree_nrays, f, lfree_rays, beta, q, + map->ray_scale[i], &lp, &value); + abort_if(rval, "GREEDY_ND_pi failed"); + } + else + { + rval = GREEDY_ND_psi(nrows, lfree_nrays, f, lfree_rays, beta, q, + map->ray_scale[i], &lp, &value); + abort_if(rval, "GREEDY_ND_psi failed"); + } + + log_verbose(" psi[%4d] = %20.12lf %d\n", map->indices[i], value); + + value *= 1.0001; + value = DOUBLE_max(value, 0.0001); + + cut->indices[i] = map->indices[i]; + cut->pi[i] = -value; + } + + cut->pi_zero = -1.0; + +CLEANUP: + LP_free(&lp); + return rval; +} + +static int select_rays(const struct RayMap map, + const struct Tableau *tableau, + struct MultiRowModel *model) +{ + int rval = 0; + int nrows = tableau->nrows; + + for (double norm_cutoff = 0.00; norm_cutoff <= 5.0; norm_cutoff += 0.1) + { + model->nrays = 0; + + for (int i = 0; i < map.nrays; i++) + { + int keep = 1; + double *r = &map.rays[tableau->nrows * i]; + + for (int j = 0; j < (model->nrays); j++) + { + double *q = &map.rays[tableau->nrows * j]; + double norm = 0; + + for (int k = 0; k < nrows; k++) + norm += fabs(r[k] - q[k]); + + if (norm <= norm_cutoff) + { + keep = 0; + break; + } + } + + if (keep) + { + memcpy(&model->rays[nrows * (model->nrays)], r, + nrows * sizeof(double)); + model->nrays++; + } + } + + log_debug(" norm_cutoff=%8.2lf nrays=%8d\n", norm_cutoff, + model->nrays); + + if (model->nrays < MAX_N_RAYS) break; + } + +CLEANUP: + return rval; +} + +static int +append_extra_rays(const struct Tableau *tableau, double *rays, int *nrays) +{ + int rval = 0; + int nrows = tableau->nrows; + int n_extra_rays = 0; + double extra_rays[100]; + + abort_if(nrows > 3, "not implemented"); + + if (nrows == 2) + { + extra_rays[0] = 0.0001; + extra_rays[1] = 0.0000; + + extra_rays[2] = -0.0001; + extra_rays[3] = 0.0001; + + extra_rays[4] = -0.0001; + extra_rays[5] = -0.0001; + + n_extra_rays = 3; + } + + if (nrows == 3) + { + extra_rays[0] = 0.0000; + extra_rays[1] = 0.0000; + extra_rays[2] = 0.0001; + + extra_rays[3] = 0.0001; + extra_rays[4] = 0.0000; + extra_rays[5] = -0.0001; + + extra_rays[6] = -0.0001; + extra_rays[7] = 0.0000; + extra_rays[8] = -0.0001; + + extra_rays[9] = -0.0001; + extra_rays[10] = -0.0001; + extra_rays[11] = -0.0001; + + n_extra_rays = 4; + } + + for (int i = 0; i < n_extra_rays; i++) + { + double *r = &extra_rays[nrows * i]; + + double scale; + int found, index; + rval = CG_find_ray(nrows, rays, *nrays, r, &found, &scale, &index); + abort_if(rval, "CG_find_ray failed"); + + if (!found) + { + memcpy(&rays[nrows * (*nrays)], r, nrows * + sizeof(double)); + (*nrays)++; + } + } + +CLEANUP: + return rval; +} + +static int write_sage_file(int nrows, + int nrays, + const double *f, + const double *rays, + const double *beta, + const char *filename) +{ + int rval = 0; + + FILE *fsage = fopen(filename, "w"); + abort_iff(!fsage, "could not open %s", filename); + + fprintf(fsage, "f=vector(["); + for (int i = 0; i < nrows; i++) + fprintf(fsage, "%.20lf,", f[i]); + fprintf(fsage, "])\n"); + + fprintf(fsage, "R=matrix([\n"); + for (int i = 0; i < nrays; i++) + { + fprintf(fsage, " ["); + for (int j = 0; j < nrows; j++) + { + fprintf(fsage, "%.20lf,", rays[i * nrows + j]); + } + fprintf(fsage, "],\n"); + } + fprintf(fsage, "])\n"); + + fprintf(fsage, "pi=vector([\n"); + for (int k = 0; k < nrays; k++) + fprintf(fsage, " %.12lf,\n", 1 / beta[k]); + fprintf(fsage, "])\n"); + +CLEANUP: + fclose(fsage); + return rval; +} + +static int dump_cut(const struct Tableau *tableau, + const double *rays, + int nrays, + const double *f, + const double *beta) +{ + int rval = 0; + + char filename[100]; + sprintf(filename, "cut-%03d.sage", DUMP_CUT_N++); + + time_printf("Writing %s...\n", filename); + rval = write_sage_file(tableau->nrows, nrays, f, rays, beta, filename); + abort_if(rval, "write_sage_file failed"); + +CLEANUP: + return rval; +} + +#ifndef TEST_SOURCE + +int INFINITY_generate_cut(struct Tableau *tableau, struct Row *cut) +{ + int rval = 0; + int max_nrays = 0; + int nrows = tableau->nrows; + double *beta = 0; + + for (int i = 0; i < tableau->nrows; i++) + max_nrays += tableau->rows[i]->nz; + + struct MultiRowModel model; + rval = CG_init_model(&model, nrows, max_nrays + 100); + abort_if(rval, "CG_init_model failed"); + + rval = CG_extract_f_from_tableau(tableau, model.f); + abort_if(rval, "CG_extract_f_from_tableau failed"); + + struct RayMap map; + rval = CG_init_ray_map(&map, max_nrays, nrows); + abort_if(rval, "CG_init_ray_map failed"); + + rval = CG_extract_rays_from_tableau(tableau, &map); + abort_if(rval, "CG_extract_rays_from_rows failed"); + + rval = select_rays(map, tableau, &model); + abort_if(rval, "select_rays failed"); + + if (ENABLE_LIFTING) + { + rval = append_extra_rays(tableau, model.rays, &model.nrays); + abort_if(rval, "append_extra_rays failed"); + } + + beta = (double *) malloc(model.nrays * sizeof(double)); + abort_if(!beta, "could not allocate beta"); + + if (model.nrays < 3) + { + rval = ERR_NO_CUT; + cut->pi = 0; + cut->indices = 0; + goto CLEANUP; + } + + log_debug("Selected %d rays\n", nrays); + if_verbose_level print_rays(tableau, &map, model.f, map.rays, map.nrays); + + log_verbose("Computing lattice-free set...\n"); + if (nrows == 2) + { + rval = sort_rays_angle(model.rays, model.nrays, beta); + abort_if(rval, "sort_rays_angle failed"); + + rval = INFINITY_2D_generate_cut(&model, beta); + if (rval) + { + rval = ERR_NO_CUT; + goto CLEANUP; + } + } + else + { + rval = GREEDY_ND_generate_cut(nrows, model.nrays, model.f, model.rays, + beta); + if (rval) + { + rval = ERR_NO_CUT; + goto CLEANUP; + } + } + + if (SHOULD_DUMP_CUTS) + { + rval = dump_cut(tableau, model.rays, model.nrays, model.f, beta); + abort_if(rval, "dump_cut failed"); + } + + rval = create_cut(tableau, &map, model.f, model.nrays, model.rays, beta, + cut); + abort_if(rval, "create_cut failed"); + + if_verbose_level print_cut(cut); + +CLEANUP: + if (beta) free(beta); + CG_free_model(&model); + CG_free_ray_map(&map); + return rval; +} + +#endif // TEST_SOURCE diff --git a/infinity/library/tests/greedy-nd-test.cpp b/infinity/library/tests/greedy-nd-test.cpp index 41b8fa2..fea9360 100644 --- a/infinity/library/tests/greedy-nd-test.cpp +++ b/infinity/library/tests/greedy-nd-test.cpp @@ -16,6 +16,8 @@ #include +#define TEST_SOURCE + extern "C" { #include #include @@ -26,7 +28,7 @@ extern "C" { int ENABLE_LIFTING = 0; int MIN_N_ROWS = 2; int MAX_N_ROWS = 2; -int DUMP_CUT = 0; +int SHOULD_DUMP_CUTS = 0; int DUMP_CUT_N = 0; TEST(GreedyNDTest, find_violated_cone_test) diff --git a/infinity/library/tests/greedy-2d-test.cpp b/infinity/library/tests/infinity-2d-test.cpp similarity index 89% rename from infinity/library/tests/greedy-2d-test.cpp rename to infinity/library/tests/infinity-2d-test.cpp index 5469eed..3609f59 100644 --- a/infinity/library/tests/greedy-2d-test.cpp +++ b/infinity/library/tests/infinity-2d-test.cpp @@ -21,22 +21,30 @@ using namespace std; extern "C" { #include -#include -#include "../src/greedy-2d.c" +#include +#include "../src/infinity-2d.c" } #define BOUNDS_EPSILON 0.01 -TEST(Greedy2DTest, test_generate_cut_1) +TEST(Infinity2DTest, test_generate_cut_1) { int rval = 0; double bounds[100]; + double f[] = {1 / 4.0, 3 / 4.0}; - double rays[] = {-2 / 5.0, 5 / 7.0, 0.0, 1.0, 1.0, 1.0, 4 / 5.0, -2 / 3.0, - -1.0, 0.0}; + double rays[] = { + -2 / 5.0, 5 / 7.0, + 0.0, 1.0, + 1.0, 1.0, + 4 / 5.0, -2 / 3.0, + -1.0, 0.0 + }; - rval = GREEDY_2D_generate_cut(rays, 5, f, bounds); - abort_if(rval, "GREEDY_2D_generate_cut failed"); + const struct MultiRowModel model = {f , rays, 5, 2}; + + rval = INFINITY_2D_generate_cut(&model, bounds); + abort_if(rval, "INFINITY_2D_generate_cut failed"); EXPECT_NEAR(23 / 50.0, bounds[0], BOUNDS_EPSILON); EXPECT_NEAR(23 / 42.0, bounds[1], BOUNDS_EPSILON); @@ -48,7 +56,7 @@ TEST(Greedy2DTest, test_generate_cut_1) if (rval) FAIL(); } -TEST(Greedy2DTest, test_generate_cut_2) +TEST(Infinity2DTest, test_generate_cut_2) { int rval = 0; double bounds[100]; @@ -61,8 +69,10 @@ TEST(Greedy2DTest, test_generate_cut_2) 1.0, -1.0 }; - rval = GREEDY_2D_generate_cut(rays, 5, f, bounds); - abort_if(rval, "GREEDY_2D_generate_cut failed"); + const struct MultiRowModel model = {f , rays, 5, 2}; + + rval = INFINITY_2D_generate_cut(&model, bounds); + abort_if(rval, "INFINITY_2D_generate_cut failed"); EXPECT_NEAR(0.5, bounds[0], BOUNDS_EPSILON); EXPECT_NEAR(0.5, bounds[1], BOUNDS_EPSILON); @@ -74,15 +84,17 @@ TEST(Greedy2DTest, test_generate_cut_2) if (rval) FAIL(); } -TEST(Greedy2DTest, test_generate_cut_3) +TEST(Infinity2DTest, test_generate_cut_3) { int rval = 0; double bounds[100]; double f[] = {5 / 22.0, 0.0}; double rays[] = {-1 / 22.0, 0.0, 0.0, 1 / 18.0, 1 / 22.0, 0.0}; - rval = GREEDY_2D_generate_cut(rays, 3, f, bounds); - abort_if(rval, "GREEDY_2D_generate_cut failed"); + const struct MultiRowModel model = {f , rays, 3, 2}; + + rval = INFINITY_2D_generate_cut(&model, bounds); + abort_if(rval, "INFINITY_2D_generate_cut failed"); EXPECT_NEAR(5.0, bounds[0], BOUNDS_EPSILON); EXPECT_NEAR(17.0, bounds[2], BOUNDS_EPSILON); @@ -92,7 +104,7 @@ TEST(Greedy2DTest, test_generate_cut_3) if (rval) FAIL(); } -TEST(Greedy2DTest, scale_to_chull_test) +TEST(Infinity2DTest, scale_to_chull_test) { int rval = 0; @@ -140,7 +152,7 @@ CLEANUP: if(rval) FAIL(); } -TEST(Greedy2DTest, scale_to_chull_test_2) +TEST(Infinity2DTest, scale_to_chull_test_2) { int rval = 0; @@ -173,7 +185,7 @@ TEST(Greedy2DTest, scale_to_chull_test_2) if(rval) FAIL(); } -TEST(Greedy2DTest, find_containing_cone_test) +TEST(Infinity2DTest, find_containing_cone_test) { int rval = 0; @@ -213,7 +225,7 @@ TEST(Greedy2DTest, find_containing_cone_test) if(rval) FAIL(); } -TEST(Greedy2DTest, find_containing_cone_test_2) +TEST(Infinity2DTest, find_containing_cone_test_2) { int rval = 0; @@ -232,7 +244,7 @@ TEST(Greedy2DTest, find_containing_cone_test_2) if(rval) FAIL(); } -TEST(Greedy2DTest, find_containing_cone_test_3) +TEST(Infinity2DTest, find_containing_cone_test_3) { int rval = 0; @@ -257,7 +269,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) if(rval) FAIL(); } -//TEST(Greedy2DTest, test_generate_cut_4) +//TEST(Infinity2DTest, test_generate_cut_4) //{ // int rval = 0; // double bounds[100]; @@ -265,8 +277,8 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // double rays[] = {0, -1 / 38.0, -1 / 22.0, -1 / 38.0, 0, 1 / 38.0, -1 / 22.0, // 0, 1 / 22.0, 0, 1 / 22.0, 1 / 38.0}; // -// rval = GREEDY_2D_generate_cut(rays, 6, f, bounds); -// abort_if(rval, "GREEDY_2D_generate_cut failed"); +// rval = INFINITY_2D_generate_cut(rays, 6, f, bounds); +// abort_if(rval, "INFINITY_2D_generate_cut failed"); // // EXPECT_NEAR(20.0, bounds[0], BOUNDS_EPSILON); // EXPECT_NEAR(20.0, bounds[1], BOUNDS_EPSILON); @@ -279,7 +291,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // if (rval) FAIL(); //} // -//TEST(Greedy2DTest, test_generate_cut_5) +//TEST(Infinity2DTest, test_generate_cut_5) //{ // int rval = 0; // double bounds[100]; @@ -291,8 +303,8 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // 0.04545454545454545581, 0.00000000000000000000, // 0.04545454545454545581, 0.02631578947368420907}; // -// rval = GREEDY_2D_generate_cut(rays, 6, f, bounds); -// abort_if(rval, "GREEDY_2D_generate_cut failed"); +// rval = INFINITY_2D_generate_cut(rays, 6, f, bounds); +// abort_if(rval, "INFINITY_2D_generate_cut failed"); // // EXPECT_NEAR(20.0, bounds[0], BOUNDS_EPSILON); // EXPECT_NEAR(20.0, bounds[1], BOUNDS_EPSILON); @@ -307,7 +319,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // // // -//TEST(Greedy2DTest, get_peak_ray_test_1) +//TEST(Infinity2DTest, get_peak_ray_test_1) //{ // int rval = 0; // @@ -340,7 +352,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // if(rval) FAIL(); //} // -//TEST(Greedy2DTest, get_peak_ray_test_2) +//TEST(Infinity2DTest, get_peak_ray_test_2) //{ // int rval = 0; // @@ -383,7 +395,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // if(rval) FAIL(); //} // -//TEST(Greedy2DTest, get_peak_ray_test_3) +//TEST(Infinity2DTest, get_peak_ray_test_3) //{ // int rval = 0; // @@ -401,7 +413,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // if(rval) FAIL(); //} // -//TEST(Greedy2DTest, nearest_lattice_point_test_1) +//TEST(Infinity2DTest, nearest_lattice_point_test_1) //{ // int rval = 0; // @@ -440,7 +452,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // if(rval) FAIL(); //} // -//TEST(Greedy2DTest, nearest_lattice_point_test_2) +//TEST(Infinity2DTest, nearest_lattice_point_test_2) //{ // int rval = 0; // @@ -460,7 +472,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // if(rval) FAIL(); //} // -//TEST(Greedy2DTest, find_normal_test) +//TEST(Infinity2DTest, find_normal_test) //{ // int rval = 0; // @@ -496,7 +508,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // if(rval) FAIL(); //} // -//TEST(Greedy2DTest, check_rays_parallel) +//TEST(Infinity2DTest, check_rays_parallel) //{ // double r1[] = {1.0, 2.0}; // double r2[] = {2.0, 4.0}; @@ -521,7 +533,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // EXPECT_FALSE(match); //} // -//TEST(Greedy2DTest, find_ray) +//TEST(Infinity2DTest, find_ray) //{ // double rays[] = {1.0, 2.0, -1.0, 0.0, -5.0, -5.0}; // int nrays = 3; @@ -557,7 +569,7 @@ TEST(Greedy2DTest, find_containing_cone_test_3) // EXPECT_FALSE(found); //} // -//TEST(Greedy2DTest, extract_rays_from_two_sparse_rows_test) +//TEST(Infinity2DTest, extract_rays_from_two_sparse_rows_test) //{ // int rval = 0; // diff --git a/infinity/library/tests/infinity-test.cpp b/infinity/library/tests/infinity-test.cpp new file mode 100644 index 0000000..812009a --- /dev/null +++ b/infinity/library/tests/infinity-test.cpp @@ -0,0 +1,91 @@ +/* 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 . + */ + +#include + +#define TEST_SOURCE + +extern "C" { +#include "../src/infinity.c" +} + +TEST(InfinityTest, cmp_ray_angle_test) +{ + double r0[] = { 1.0, 0.0 }; + double r1[] = { 2.0, 0.0 }; + double r2[] = { 1.0, 1.0 }; + double r3[] = { -1.0, 0.0 }; + double r4[] = { 1.0, -1.0 }; + + SortPair sp0 = { 0, &r0 }; + SortPair sp1 = { 1, &r1 }; + SortPair sp2 = { 2, &r2 }; + SortPair sp3 = { 3, &r3 }; + SortPair sp4 = { 4, &r4 }; + + EXPECT_EQ(_qsort_cmp_rays_angle(&sp0, &sp1), 0); + EXPECT_EQ(_qsort_cmp_rays_angle(&sp0, &sp2), 1); + EXPECT_EQ(_qsort_cmp_rays_angle(&sp0, &sp3), 1); + EXPECT_EQ(_qsort_cmp_rays_angle(&sp0, &sp4), -1); + + EXPECT_EQ(_qsort_cmp_rays_angle(&sp2, &sp0), -1); + EXPECT_EQ(_qsort_cmp_rays_angle(&sp2, &sp1), -1); + EXPECT_EQ(_qsort_cmp_rays_angle(&sp2, &sp3), 1); + EXPECT_EQ(_qsort_cmp_rays_angle(&sp2, &sp4), -1); + + EXPECT_EQ(_qsort_cmp_rays_angle(&sp3, &sp0), -1); +} + +TEST(InfinityTest, sort_rays_angle_test) +{ + int rval = 0; + + int n_rays = 5; + double beta[] = {0, 1, 2, 3, 4}; + SortPair sp0, sp1, sp2, sp3, sp4; + + double rays[] = { + 1.0, 1.0, + 1.0, 0.0, + 1.0, -1.0, + -1.0, 0.0, + 2.0, 0.0, + }; + + rval = sort_rays_angle(rays, n_rays, beta); + abort_if(rval, "sort_rays_angle failed"); + + sp0 = { 0, &rays[0] }; + sp1 = { 1, &rays[2] }; + sp2 = { 2, &rays[4] }; + sp3 = { 3, &rays[6] }; + sp4 = { 4, &rays[8] }; + + EXPECT_LE(_qsort_cmp_rays_angle(&sp0, &sp1), 0); + EXPECT_LE(_qsort_cmp_rays_angle(&sp1, &sp2), 0); + EXPECT_LE(_qsort_cmp_rays_angle(&sp2, &sp3), 0); + EXPECT_LE(_qsort_cmp_rays_angle(&sp3, &sp4), 0); + + EXPECT_EQ(beta[0], 3); + EXPECT_EQ(beta[1], 0); + EXPECT_TRUE(beta[2] == 1 || beta[2] == 4); + EXPECT_TRUE(beta[3] == 1 || beta[3] == 4); + EXPECT_TRUE(beta[2] != beta[3]); + EXPECT_EQ(beta[4], 2); + +CLEANUP: + if (rval) FAIL(); +} diff --git a/multirow/include/multirow/cg.h b/multirow/include/multirow/cg.h index db2a07a..e373077 100644 --- a/multirow/include/multirow/cg.h +++ b/multirow/include/multirow/cg.h @@ -38,42 +38,54 @@ struct CG double *current_solution; }; +struct Tableau +{ + int nrows; + struct Row **rows; + char *column_types; +}; + +struct RayMap +{ + int nrays; + double *rays; + int *variable_to_ray; + double *ray_scale; + int *indices; + int nvars; +}; + +struct MultiRowModel +{ + double *f; + double *rays; + int nrays; + int nrows; +}; + typedef int (*SingleRowGeneratorCallback)(const struct Row *row, char *column_types, struct Row *cut); -typedef int (*MultirowGeneratorCallback)(int nrows, - struct Row **rows, - char *column_types, +typedef int (*MultiRowGeneratorCallback)(const struct Tableau *tableau, struct Row *cut); -int CG_init(struct LP *lp, - char *column_types, - struct CG *cg); +int CG_init(struct LP *lp, char *column_types, struct CG *cg); void CG_free(struct CG *cg); -int CG_add_single_row_cuts(struct CG *cg, - SingleRowGeneratorCallback generate); +int CG_add_single_row_cuts(struct CG *cg, SingleRowGeneratorCallback generate); int CG_add_multirow_cuts(struct CG *cg, int nrows, - MultirowGeneratorCallback generate); + MultiRowGeneratorCallback generate); -int CG_set_integral_solution(struct CG *cg, - double *valid_solution); +int CG_set_integral_solution(struct CG *cg, double *valid_solution); -int CG_set_basic_solution(struct CG *cg, - double *basic_solution); +int CG_set_basic_solution(struct CG *cg, double *basic_solution); -int CG_extract_rays_from_rows(int nrows, - struct Row **rows, - double *rays, - int *nrays, - int *variable_to_ray, - double *ray_scale, - int *indices, - int *nz); +int CG_extract_rays_from_tableau(const struct Tableau *tableau, + struct RayMap *map); int CG_boost_variable(int var, double factor, @@ -92,4 +104,14 @@ int CG_find_ray(int dim, double *scale, int *index); +int CG_init_ray_map(struct RayMap *map, int max_nrays, int nrows); + +void CG_free_ray_map(struct RayMap *map); + +int CG_extract_f_from_tableau(const struct Tableau *tableau, double *f); + +int CG_init_model(struct MultiRowModel *model, int nrows, int max_nrays); + +void CG_free_model(struct MultiRowModel *model); + #endif //MULTIROW_CG_H diff --git a/multirow/include/multirow/params.h b/multirow/include/multirow/params.h index 7a48be7..55be1bb 100644 --- a/multirow/include/multirow/params.h +++ b/multirow/include/multirow/params.h @@ -77,7 +77,7 @@ extern int BOOST_VAR; extern double BOOST_FACTOR; -extern int DUMP_CUT; +extern int SHOULD_DUMP_CUTS; extern int DUMP_CUT_N; extern int ENABLE_LIFTING; diff --git a/multirow/src/cg.c b/multirow/src/cg.c index 3f3bc1f..6911c56 100644 --- a/multirow/src/cg.c +++ b/multirow/src/cg.c @@ -31,27 +31,27 @@ static int select_rows(struct CG *cg, int *row_selected) long histogram[100] = {0}; - for(int i = 0; i < cg->nrows; i++) + for (int i = 0; i < cg->nrows; i++) { row_selected[i] = 1; struct Row *r1 = cg->tableau_rows[i]; - if(r1->head < 0 || cg->column_types[r1->head] != MILP_INTEGER) + if (r1->head < 0 || cg->column_types[r1->head] != MILP_INTEGER) { row_selected[i] = 0; continue; } double dist = fabs(0.5 - frac(r1->pi_zero)); - for(int k = 0; k < 100; k++) - if(dist <= (k / 100.0)) - histogram[k]++; + for (int k = 0; k < 100; k++) + if (dist <= (k / 100.0)) + histogram[k]++; } double dist_cutoff = 0.5; - for(int k = 0; k < 100; k++) + for (int k = 0; k < 100; k++) { - if(histogram[k] > MAX_SELECTED_ROWS) + if (histogram[k] > MAX_SELECTED_ROWS) { dist_cutoff = k / 100.0; break; @@ -59,14 +59,14 @@ static int select_rows(struct CG *cg, int *row_selected) } int selected_count = 0; - for(int i = 0; i < cg->nrows; i++) + for (int i = 0; i < cg->nrows; i++) { - if(!row_selected[i]) continue; + if (!row_selected[i]) continue; struct Row *r1 = cg->tableau_rows[i]; double dist = fabs(0.5 - frac(r1->pi_zero)); - if(dist > dist_cutoff || selected_count >= MAX_SELECTED_ROWS) + if (dist > dist_cutoff || selected_count >= MAX_SELECTED_ROWS) { row_selected[i] = 0; continue; @@ -79,11 +79,8 @@ CLEANUP: return rval; } -static int next_combination(int n, - int k, - int *array, - int inc_index, - int *finished) +static int +next_combination(int n, int k, int *array, int inc_index, int *finished) { int i; int rval = 0; @@ -91,13 +88,13 @@ static int next_combination(int n, array[inc_index]++; *finished = 0; - for(i = inc_index; i < k; i++) + for (i = inc_index; i < k; i++) { - if(array[i] < n - i) + if (array[i] < n - i) break; - if(i+1 < k) - array[i+1]++; + if (i + 1 < k) + array[i + 1]++; else { *finished = 1; @@ -105,26 +102,22 @@ static int next_combination(int n, } } - for(int j = 1; j <= i; j++) + for (int j = 1; j <= i; j++) array[i - j] = array[i] + j; return 0; } - -static int ray_norm(int dim, - const double *ray, - double *norm) +static int ray_norm(int dim, const double *ray, double *norm) { *norm = 0; - for(int i = 0; i < dim; i++) + for (int i = 0; i < dim; i++) *norm += fabs(ray[i]); return 0; } - /* * Checks whether two rays are parallel. Also returns the scaling factor, * in case they are. @@ -157,16 +150,17 @@ static int check_rays_parallel(int dim, } else { - for(int i = 0; i < dim; i++) + for (int i = 0; i < dim; i++) { - if(DOUBLE_sgn(r1[i]) != DOUBLE_sgn(r2[i]) || + if (DOUBLE_sgn(r1[i]) != DOUBLE_sgn(r2[i]) || DOUBLE_neq(r1[i] * r2_norm, r2[i] * r1_norm)) { *match = 0; break; } - log_verbose(" %.12lf equals %.12lf\n", r1[i] * r2_norm, r2[i] * r1_norm); + log_verbose(" %.12lf equals %.12lf\n", r1[i] * r2_norm, + r2[i] * r1_norm); } } @@ -176,9 +170,7 @@ CLEANUP: return rval; } - -static void print_row(const struct CG *cg, - const struct Row *row) +static void print_row(const struct CG *cg, const struct Row *row) { double *x = cg->integral_solution; double *y = cg->basic_solution; @@ -195,10 +187,8 @@ static void print_row(const struct CG *cg, time_printf(" <= %20.6lf [%d]\n", row->pi_zero, row->head); } - -static int evaluate_row_pair(const struct Row *row1, - const struct Row *row2, - double *score) +static int +evaluate_row_pair(const struct Row *row1, const struct Row *row2, double *score) { int rval = 0; @@ -217,7 +207,7 @@ static int evaluate_row_pair(const struct Row *row1, if (i2 < row2->nz) idx2 = row2->indices[i2]; - if(idx1 == idx2) hit++; + if (idx1 == idx2) hit++; int idx_min = min(idx1, idx2); @@ -230,14 +220,12 @@ static int evaluate_row_pair(const struct Row *row1, *score = (hit * 2.0) / (row1->nz + row2->nz); - CLEANUP: +CLEANUP: return rval; } -static double replace_x(const double *pi, - const int *indices, - int nz, - const double *x) +static double +replace_x(const double *pi, const int *indices, int nz, const double *x) { double lhs = 0; @@ -247,7 +235,7 @@ static double replace_x(const double *pi, int idx = indices[i]; if (!DOUBLE_iszero(pii) && !DOUBLE_iszero(x[idx])) - log_verbose(" %12.8lf * %12.8lf (x%d)\n", pii, x[idx], idx); + log_verbose(" %12.8lf * %12.8lf (x%d)\n", pii, x[idx], idx); lhs += pii * x[idx]; } @@ -255,10 +243,7 @@ static double replace_x(const double *pi, return lhs; } - -static int copy_solution(struct CG *cg, - double *from, - double **to) +static int copy_solution(struct CG *cg, double *from, double **to) { int rval = 0; int ncols = LP_get_num_cols(cg->lp); @@ -275,9 +260,7 @@ CLEANUP: return rval; } - -static int check_cut(struct CG *cg, - struct Row *cut) +static int check_cut(struct CG *cg, struct Row *cut) { int rval = 0; @@ -286,7 +269,7 @@ static int check_cut(struct CG *cg, time_printf("Checking cut:\n"); for (int i = 0; i < cut->nz; i++) { - if(DOUBLE_iszero(cg->integral_solution[cut->indices[i]])) continue; + if (DOUBLE_iszero(cg->integral_solution[cut->indices[i]])) continue; time_printf(" %12.8lf x%d", cut->pi[i], cut->indices[i]); if (cg->integral_solution) printf(" (=%12.8lf)", cg->integral_solution[cut->indices[i]]); @@ -300,11 +283,11 @@ static int check_cut(struct CG *cg, log_verbose("Basic solution check:\n"); double lhs = replace_x(cut->pi, cut->indices, cut->nz, - cg->basic_solution); + cg->basic_solution); log_verbose(" %.8lf > %.8lf\n", lhs, cut->pi_zero); abort_iff(!DOUBLE_geq(lhs, cut->pi_zero), "Cut fails to cut " - "basic solution: %12.8lf < %12.8lf", lhs, cut->pi_zero); + "basic solution: %12.8lf < %12.8lf", lhs, cut->pi_zero); } if (cg->integral_solution) @@ -312,7 +295,7 @@ static int check_cut(struct CG *cg, log_verbose("Integral solution check:\n"); double lhs = replace_x(cut->pi, cut->indices, cut->nz, - cg->integral_solution); + cg->integral_solution); log_verbose(" %.8lf <= %.8lf\n", lhs, cut->pi_zero); abort_iff(!DOUBLE_leq(lhs, cut->pi_zero), "Cut cuts off known integral " @@ -323,10 +306,7 @@ CLEANUP: return rval; } - -static int add_cut(struct CG *cg, - struct Row *cut, - int *ignored) +static int add_cut(struct CG *cg, struct Row *cut, int *ignored) { int rval; double *x = 0; @@ -338,13 +318,12 @@ static int add_cut(struct CG *cg, rval = check_cut(cg, cut); abort_if(rval, "check_cut failed"); - lhs = replace_x(cut->pi, cut->indices, cut->nz, - cg->current_solution); + lhs = replace_x(cut->pi, cut->indices, cut->nz, cg->current_solution); *ignored = 0; STATS_increment_generated_cuts(); - if(DOUBLE_leq(lhs, cut->pi_zero)) + if (DOUBLE_leq(lhs, cut->pi_zero)) { log_verbose("Ignoring cut (%12.8lf <= %12.8lf)\n", lhs, cut->pi_zero); *ignored = 1; @@ -363,12 +342,12 @@ static int add_cut(struct CG *cg, rval = LP_get_obj_val(cg->lp, &obj); abort_if(rval, "LP_get_obj_val failed"); - if(DOUBLE_neq(obj, cg->last_obj_value)) + if (DOUBLE_neq(obj, cg->last_obj_value)) log_info(" opt = %lf\n", obj); cg->last_obj_value = obj; - x = (double*) malloc(cg->ncols * sizeof(double)); + x = (double *) malloc(cg->ncols * sizeof(double)); abort_if(!x, "could not allocate x"); rval = LP_get_x(cg->lp, x); @@ -381,7 +360,7 @@ static int add_cut(struct CG *cg, } CLEANUP: - if(x) free(x); + if (x) free(x); return rval; } @@ -419,43 +398,40 @@ CLEANUP: return 0; } -int CG_extract_rays_from_rows(int nrows, - struct Row **rows, - double *rays, - int *nrays, - int *variable_to_ray, - double *ray_scale, - int *indices, - int *nz) +int CG_extract_rays_from_tableau(const struct Tableau *tableau, + struct RayMap *map) { int rval = 0; - (*nz) = 0; - (*nrays) = 0; + int nrows = tableau->nrows; + double *rays = map->rays; + struct Row **rows = tableau->rows; + + map->nvars = 0; + map->nrays = 0; int *i = 0; int *idx = 0; - - i = (int*) malloc(nrows * sizeof(int)); - idx = (int*) malloc(nrows * sizeof(int)); + i = (int *) malloc(nrows * sizeof(int)); + idx = (int *) malloc(nrows * sizeof(int)); abort_if(!i, "could not allocate i"); abort_if(!idx, "could not allocate idx"); - for(int j = 0; j < nrows; j++) + for (int j = 0; j < nrows; j++) i[j] = 0; while (1) { - double *r = &rays[nrows * (*nrays)]; + double *r = &rays[nrows * (map->nrays)]; int idx_min = INT_MAX; - for(int j = 0; j < nrows; j++) + for (int j = 0; j < nrows; j++) { r[j] = 0.0; - if(i[j] < rows[j]->nz) + if (i[j] < rows[j]->nz) idx[j] = rows[j]->indices[i[j]]; else idx[j] = INT_MAX; @@ -463,18 +439,18 @@ int CG_extract_rays_from_rows(int nrows, idx_min = min(idx_min, idx[j]); } - if(idx_min == INT_MAX) + if (idx_min == INT_MAX) break; - for(int j = 0; j < nrows; j++) + for (int j = 0; j < nrows; j++) { - if(idx[j] > idx_min) continue; + if (idx[j] > idx_min) continue; r[j] = -rows[j]->pi[i[j]]; i[j]++; } - for(int j = 0; j < nrows; j++) - if(idx_min == rows[j]->head) + for (int j = 0; j < nrows; j++) + if (idx_min == rows[j]->head) goto NEXT_RAY; int found; @@ -483,67 +459,65 @@ int CG_extract_rays_from_rows(int nrows, log_verbose(" extracted ray (%d):\n", idx_min); - for(int j=0; j < nrows; j++) - log_verbose(" r[%d] = %.12lf\n", j, r[j]); + for (int j = 0; j < nrows; j++) + log_verbose(" r[%d] = %.12lf\n", j, r[j]); - rval = CG_find_ray(nrows, rays, *nrays, r, &found, &scale, &ray_index); + rval = CG_find_ray(nrows, rays, map->nrays, r, &found, &scale, + &ray_index); abort_if(rval, "CG_find_ray failed"); if (!found) { log_verbose(" ray is new\n"); scale = 1.0; - ray_index = (*nrays)++; + ray_index = (map->nrays)++; } else { log_verbose(" ray equals:\n"); double *q = &rays[ray_index * nrows]; - for(int j=0; j < nrows; j++) - log_verbose(" r[%d] = %.12lf\n", j, q[j]); + for (int j = 0; j < nrows; j++) + log_verbose(" r[%d] = %.12lf\n", j, q[j]); } - ray_scale[(*nz)] = scale; - indices[(*nz)] = idx_min; - variable_to_ray[(*nz)] = ray_index; - (*nz)++; + map->ray_scale[map->nvars] = scale; + map->indices[map->nvars] = idx_min; + map->variable_to_ray[map->nvars] = ray_index; + map->nvars++; - NEXT_RAY: - ; + NEXT_RAY:; } - for(int j = 0; j < *nrays; j++) + for (int j = 0; j < map->nrays; j++) { double *r = &rays[nrows * j]; double max_scale = 0.0; - for(int k = 0; k < *nz; k++) + for (int k = 0; k < map->nvars; k++) { - if(variable_to_ray[k] != j) continue; - if(ray_scale[k] < max_scale) continue; - max_scale = ray_scale[k]; + if (map->variable_to_ray[k] != j) continue; + if (map->ray_scale[k] < max_scale) continue; + max_scale = map->ray_scale[k]; } abort_if(max_scale == 0.0, "max_scale is zero"); - for(int k = 0; k < *nz; k++) - if(variable_to_ray[k] == j) ray_scale[k] /= max_scale; + for (int k = 0; k < map->nvars; k++) + if (map->variable_to_ray[k] == j) + map->ray_scale[k] /= max_scale; - for(int k = 0; k < nrows; k++) + for (int k = 0; k < nrows; k++) r[k] *= max_scale; } CLEANUP: - if(idx) free(idx); - if(i) free(i); + if (idx) free(idx); + if (i) free(i); return rval; } - -int CG_init(struct LP *lp, - char *column_types, - struct CG *cg) +int CG_init(struct LP *lp, char *column_types, struct CG *cg) { int rval = 0; @@ -581,7 +555,7 @@ int CG_init(struct LP *lp, abort_if(!cg->tableau_rows, "could not allocate cg->tableau_rows"); rval = LP_get_tableau(lp, cg->tableau_rows, cg->cstat, cg->rstat, cg->ub, - cg->lb); + cg->lb); abort_if(rval, "LP_get_tableau failed"); cg_initial_time = get_user_time(); @@ -590,7 +564,6 @@ CLEANUP: return rval; } - void CG_free(struct CG *cg) { if (!cg) return; @@ -613,9 +586,7 @@ void CG_free(struct CG *cg) free(cg); } - -int CG_add_single_row_cuts(struct CG *cg, - SingleRowGeneratorCallback generate) +int CG_add_single_row_cuts(struct CG *cg, SingleRowGeneratorCallback generate) { int rval = 0; @@ -663,80 +634,79 @@ int estimate_multirow_cut_count(struct CG *cg, *total_count = 0; - row_indices = (int*) malloc(nrows * sizeof(int)); + row_indices = (int *) malloc(nrows * sizeof(int)); abort_if(!row_indices, "could not allocate row_indices"); - for(int i = 0; i < nrows; i++) + for (int i = 0; i < nrows; i++) row_indices[i] = nrows - i - 1; - for(int i = 0; i < cg->nrows * cg->nrows; i++) + for (int i = 0; i < cg->nrows * cg->nrows; i++) row_affinity[i] = -1; - + do { int inc_index = 0; int is_rhs_integer = 1; int valid_combination = 1; - for(int i = 0; i < nrows; i++) + for (int i = 0; i < nrows; i++) { struct Row *r = cg->tableau_rows[row_indices[i]]; - if(!row_selected[row_indices[i]]) + if (!row_selected[row_indices[i]]) { valid_combination = 0; inc_index = i; } double df = fabs(frac(r->pi_zero) - 0.5); - if(df < INTEGRALITY_THRESHOLD) is_rhs_integer = 0; + if (df < INTEGRALITY_THRESHOLD) is_rhs_integer = 0; } - if(is_rhs_integer) valid_combination = 0; + if (is_rhs_integer) valid_combination = 0; - for(int i = 0; valid_combination && i < nrows; i++) + for (int i = 0; valid_combination && i < nrows; i++) { - for(int j = i+1; valid_combination && j < nrows; j++) + for (int j = i + 1; valid_combination && j < nrows; j++) { - int i1 = row_indices[i]; - int i2 = row_indices[j]; - if(i2 < i1) swap(i1, i2, int); - - int k = cg->nrows * i1 + i2; - - if(row_affinity[k] < 0) - { - struct Row *row1 = cg->tableau_rows[i1]; - struct Row *row2 = cg->tableau_rows[i2]; - double score; - - rval = evaluate_row_pair(row1, row2, &score); - abort_if(rval, "evaluate_row_pair failed"); - - row_affinity[k] = 1; - if(score < score_cutoff) row_affinity[k] = 0; - } - - if(!row_affinity[k]) - { - valid_combination = 0; - goto NEXT_COMBINATION; - } + int i1 = row_indices[i]; + int i2 = row_indices[j]; + if (i2 < i1) swap(i1, i2, int); + + int k = cg->nrows * i1 + i2; + + if (row_affinity[k] < 0) + { + struct Row *row1 = cg->tableau_rows[i1]; + struct Row *row2 = cg->tableau_rows[i2]; + double score; + + rval = evaluate_row_pair(row1, row2, &score); + abort_if(rval, "evaluate_row_pair failed"); + + row_affinity[k] = 1; + if (score < score_cutoff) row_affinity[k] = 0; + } + + if (!row_affinity[k]) + { + valid_combination = 0; + goto NEXT_COMBINATION; + } } } - if(valid_combination) + if (valid_combination) { (*total_count)++; } NEXT_COMBINATION: next_combination(cg->nrows, nrows, row_indices, inc_index, &finished); - } - while(!finished); + } while (!finished); CLEANUP: - if(row_indices) free(row_indices); + if (row_indices) free(row_indices); return rval; } @@ -756,10 +726,9 @@ int cut_dynamism(struct Row *cut, double *dynamism) return 0; } - int CG_add_multirow_cuts(struct CG *cg, int nrows, - MultirowGeneratorCallback generate) + MultiRowGeneratorCallback generate) { int rval = 0; int *row_indices = 0; @@ -771,10 +740,10 @@ int CG_add_multirow_cuts(struct CG *cg, long count = 0; long total_count = 0; - row_indices = (int*) malloc(nrows * sizeof(int)); + row_indices = (int *) malloc(nrows * sizeof(int)); rows = (struct Row **) malloc(nrows * sizeof(struct Row *)); - row_affinity = (int*) malloc((cg->nrows * cg->nrows) * sizeof(int)); - row_selected = (int*) malloc(cg->nrows * sizeof(int)); + row_affinity = (int *) malloc((cg->nrows * cg->nrows) * sizeof(int)); + row_selected = (int *) malloc(cg->nrows * sizeof(int)); abort_if(!row_indices, "could not allocate row_indices"); abort_if(!rows, "could not allocate rows"); @@ -785,24 +754,25 @@ int CG_add_multirow_cuts(struct CG *cg, abort_if(rval, "select_rows failed"); log_info(" Finding combinations...\n"); - for(double cutoff = 0.05; cutoff <= 1.0; cutoff += 0.05) + for (double cutoff = 0.05; cutoff <= 1.0; cutoff += 0.05) { double cg_current_time = get_user_time() - cg_initial_time; - if(cg_current_time > CG_TIMEOUT) break; + if (cg_current_time > CG_TIMEOUT) break; - rval = estimate_multirow_cut_count(cg, nrows, row_selected, row_affinity, &total_count, cutoff); + rval = estimate_multirow_cut_count(cg, nrows, row_selected, + row_affinity, &total_count, cutoff); abort_if(rval, "estimate_two_row_cuts_count failed"); log_debug(" %8d combinations [%.2lf]\n", total_count, cutoff); - if(total_count < MAX_SELECTED_COMBINATIONS) + if (total_count < MAX_SELECTED_COMBINATIONS) { log_info(" %8d combinations [%.2lf]\n", total_count, cutoff); break; } } - for(int i = 0; i < nrows; i++) + for (int i = 0; i < nrows; i++) row_indices[i] = nrows - i - 1; total_count = min(total_count, MAX_SELECTED_COMBINATIONS); @@ -813,95 +783,96 @@ int CG_add_multirow_cuts(struct CG *cg, do { double cg_current_time = get_user_time() - cg_initial_time; - if(cg_current_time > CG_TIMEOUT) break; + if (cg_current_time > CG_TIMEOUT) break; int inc_index = 0; int is_rhs_integer = 1; int valid_combination = 1; - for(int i = 0; i < nrows; i++) + for (int i = 0; i < nrows; i++) { rows[i] = cg->tableau_rows[row_indices[i]]; - if(!row_selected[row_indices[i]]) + if (!row_selected[row_indices[i]]) { valid_combination = 0; inc_index = i; } double df = fabs(frac(rows[i]->pi_zero) - 0.5); - if(df < INTEGRALITY_THRESHOLD) is_rhs_integer = 0; + if (df < INTEGRALITY_THRESHOLD) is_rhs_integer = 0; } - if(is_rhs_integer) valid_combination = 0; + if (is_rhs_integer) valid_combination = 0; - for(int i = 0; valid_combination && i < nrows; i++) + for (int i = 0; valid_combination && i < nrows; i++) { - for(int j = i+1; valid_combination && j < nrows; j++) + for (int j = i + 1; valid_combination && j < nrows; j++) { - int i1 = row_indices[i]; - int i2 = row_indices[j]; - if(i2 < i1) swap(i1, i2, int); + int i1 = row_indices[i]; + int i2 = row_indices[j]; + if (i2 < i1) swap(i1, i2, int); - int k = cg->nrows * i1 + i2; + int k = cg->nrows * i1 + i2; - abort_if(row_affinity[k] < 0, "row_affinity not computed"); - if(!row_affinity[k]) valid_combination = 0; + abort_if(row_affinity[k] < 0, "row_affinity not computed"); + if (!row_affinity[k]) valid_combination = 0; - if_verbose_level - { - struct Row *row1 = cg->tableau_rows[i1]; - struct Row *row2 = cg->tableau_rows[i2]; - double score; + if_verbose_level + { + struct Row *row1 = cg->tableau_rows[i1]; + struct Row *row2 = cg->tableau_rows[i2]; + double score; - rval = evaluate_row_pair(row1, row2, &score); - abort_if(rval, "evaluate_row_pair failed"); + rval = evaluate_row_pair(row1, row2, &score); + abort_if(rval, "evaluate_row_pair failed"); - log_verbose("%4d %4d %.2lf\n", i1, i2, score); - } + log_verbose("%4d %4d %.2lf\n", i1, i2, score); + } } } - if(valid_combination) + if (valid_combination) { count++; - if(count > MAX_SELECTED_COMBINATIONS) + if (count > MAX_SELECTED_COMBINATIONS) { - log_info(" maximum number of combinations reached. stopping.\n"); + log_info( + " maximum number of combinations reached. stopping.\n"); break; } - if_debug_level - if(ONLY_CUT > 0 && count != ONLY_CUT) + if_debug_level if (ONLY_CUT > 0 && count != ONLY_CUT) goto NEXT_COMBINATION; if_debug_level { time_printf("Generating cut %d from [ ", count); - for(int i = 0; i < nrows; i++) + for (int i = 0; i < nrows; i++) printf("%d ", row_indices[i]); printf("]...\n"); } - if(LOG_LEVEL == LOG_LEVEL_INFO) + if (LOG_LEVEL == LOG_LEVEL_INFO) { progress_print(); progress_increment(); } - + struct Row *cut = 0; - - cut = (struct Row*) malloc(sizeof(struct Row)); + + cut = (struct Row *) malloc(sizeof(struct Row)); abort_if(!cut, "could not allocate cut"); cut->pi = 0; cut->indices = 0; double initial_time = get_user_time(); + struct Tableau tableau = {nrows, rows, cg->column_types}; - rval = generate(nrows, rows, cg->column_types, cut); - if(rval == ERR_NO_CUT) + rval = generate(&tableau, cut); + if (rval == ERR_NO_CUT) { rval = 0; log_verbose("combination does not yield cut\n"); @@ -917,7 +888,7 @@ int CG_add_multirow_cuts(struct CG *cg, rval = cut_dynamism(cut, &dynamism); abort_if(rval, "cut_dynamism failed"); - if(dynamism > MAX_CUT_DYNAMISM) + if (dynamism > MAX_CUT_DYNAMISM) { log_debug("Discarding cut (dynamism=%.2lf)\n", dynamism); LP_free_row(cut); @@ -926,7 +897,7 @@ int CG_add_multirow_cuts(struct CG *cg, int ignored; rval = add_cut(cg, cut, &ignored); - if(rval) + if (rval) { log_warn("invalid cut skipped (cut %d)\n", count); rval = 0; @@ -936,39 +907,34 @@ int CG_add_multirow_cuts(struct CG *cg, } - if_debug_level if(!ignored) - { - DUMP_CUT = 1; - generate(nrows, rows, cg->column_types, cut); - DUMP_CUT = 0; - } + if_debug_level if (!ignored) + { + SHOULD_DUMP_CUTS = 1; + generate(&tableau, cut); + SHOULD_DUMP_CUTS = 0; + } LP_free_row(cut); } NEXT_COMBINATION: next_combination(cg->nrows, nrows, row_indices, inc_index, &finished); - } - while(!finished); + } while (!finished); CLEANUP: - if(row_selected) free(row_selected); - if(row_affinity) free(row_affinity); - if(row_indices) free(row_indices); - if(rows) free(rows); + if (row_selected) free(row_selected); + if (row_affinity) free(row_affinity); + if (row_indices) free(row_indices); + if (rows) free(rows); return rval; } - -int CG_set_integral_solution(struct CG *cg, - double *valid_solution) +int CG_set_integral_solution(struct CG *cg, double *valid_solution) { return copy_solution(cg, valid_solution, &cg->integral_solution); } - -int CG_set_basic_solution(struct CG *cg, - double *basic_solution) +int CG_set_basic_solution(struct CG *cg, double *basic_solution) { int rval = 0; @@ -994,28 +960,88 @@ int CG_boost_variable(int var, int rval = 0; int selected_ray = -1; - for(int j = 0; j < nz; j++) + for (int j = 0; j < nz; j++) { - if(indices[j] == var) - { - selected_ray = variable_to_ray[j]; - break; - } + if (indices[j] == var) + { + selected_ray = variable_to_ray[j]; + break; + } } - if(selected_ray < 0) goto CLEANUP; + if (selected_ray < 0) goto CLEANUP; double *r = &rays[nrows * selected_ray]; - for(int k = 0; k < nrows; k++) + for (int k = 0; k < nrows; k++) r[k] *= factor; - for(int j = 0; j < nz; j++) + for (int j = 0; j < nz; j++) { - if(variable_to_ray[j] == selected_ray) + if (variable_to_ray[j] == selected_ray) ray_scale[j] /= factor; } CLEANUP: return rval; } + +int CG_init_ray_map(struct RayMap *map, int max_nrays, int nrows) +{ + int rval = 0; + + map->variable_to_ray = (int *) malloc(max_nrays * sizeof(int)); + map->indices = (int *) malloc(max_nrays * sizeof(int)); + map->ray_scale = (double *) malloc(max_nrays * sizeof(double)); + map->rays = (double *) malloc(nrows * max_nrays * sizeof(double)); + abort_if(!map->variable_to_ray, "could not allocate variable_to_ray"); + abort_if(!map->indices, "could not allocate indices"); + abort_if(!map->ray_scale, "could not allocate ray_scale"); + abort_if(!map->rays, "could not allocate rays"); + +CLEANUP: + return rval; +} + +void CG_free_ray_map(struct RayMap *map) +{ + if(!map) return; + free(map->variable_to_ray); + free(map->indices); + free(map->ray_scale); + free(map->rays); +} + +int CG_extract_f_from_tableau(const struct Tableau *tableau, double *f) +{ + for (int i = 0; i < tableau->nrows; i++) + { + f[i] = frac(tableau->rows[i]->pi_zero); + if (DOUBLE_eq(f[i], 1.0)) f[i] = 0.0; + } + + return 0; +} + +int CG_init_model(struct MultiRowModel *model, int nrows, int max_nrays) +{ + int rval = 0; + + model->nrays = 0; + model->nrows = nrows; + model->f = (double*) malloc(nrows * sizeof(double)); + model->rays = (double*) malloc(max_nrays * sizeof(double)); + abort_if(!model->f, "could not allocate f"); + abort_if(!model->rays, "could not allocate rays"); + +CLEANUP: + return rval; +} + +void CG_free_model(struct MultiRowModel *model) +{ + if(!model) return; + free(model->rays); + free(model->f); +} + #endif // TEST_SOURCE diff --git a/multirow/tests/cg-test.cpp b/multirow/tests/cg-test.cpp index 7ffe00f..18d6524 100644 --- a/multirow/tests/cg-test.cpp +++ b/multirow/tests/cg-test.cpp @@ -151,8 +151,8 @@ TEST(CGTest, extract_rays_from_rows_test) double rays[1000]; double ray_scale[1000]; - rval = CG_extract_rays_from_rows(3, rows, rays, &nrays, variable_to_ray, - ray_scale, indices, &nz); + rval = CG_extract_rays_from_tableau(3, rows, rays, &nrays, variable_to_ray, + ray_scale, indices, &nz); abort_if(rval, "CG_extract_rays_from_rows failed"); EXPECT_EQ(nrays, 4);