diff --git a/infinity/library/CMakeLists.txt b/infinity/library/CMakeLists.txt index 0ef1fe3..2b35c4f 100644 --- a/infinity/library/CMakeLists.txt +++ b/infinity/library/CMakeLists.txt @@ -1,21 +1,19 @@ set(COMMON_SOURCES - src/infinity-2d.c - src/greedy-nd.c - src/greedy-bsearch.c - src/infinity.c - include/infinity/infinity-2d.h - include/infinity/greedy-nd.h - include/infinity/greedy-bsearch.h - include/infinity/infinity.h) + src/infinity-2d.c + src/infinity-nd.c + src/infinity.c + include/infinity/infinity-2d.h + include/infinity/infinity-nd.h + include/infinity/infinity.h) set(TEST_SOURCES - tests/infinity-2d-test.cpp - tests/greedy-nd-test.cpp - tests/infinity-test.cpp) + tests/infinity-2d-test.cpp + tests/infinity-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}) +target_include_directories(infinity_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) add_executable(infinity-test.run ${COMMON_SOURCES} ${TEST_SOURCES}) target_link_libraries(infinity-test.run gtest_main multirow_static lifting_static) diff --git a/infinity/library/include/infinity/greedy-nd.h b/infinity/library/include/infinity/greedy-nd.h deleted file mode 100644 index d37c4b8..0000000 --- a/infinity/library/include/infinity/greedy-nd.h +++ /dev/null @@ -1,74 +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 . - */ -#ifndef MULTIROW_GREEDY_ND_H -#define MULTIROW_GREEDY_ND_H - -int GREEDY_create_psi_lp(const struct ConvLFreeSet *lfree, struct LP *lp); - -int GREEDY_ND_psi(const int nrows, - const double *q, - const double q_scale, - struct LP *lp, - double *value); - -int GREEDY_ND_pi(const int nrows, - const double *q, - const double q_scale, - struct LP *lp, - double *value); - -int INFINITY_ND_generate_lfree(const struct MultiRowModel *model, - struct ConvLFreeSet *lfree); - -int GREEDY_ND_cone_bound(int nrows, - int nrays, - const double *f, - const double *rays, - const int *rx, - const double *x, - const double *beta, - double *epsilon); - -int GREEDY_ND_find_violated_cone(int nrows, - int nrays, - const double *f, - const double *rays, - const double *x, - const double *beta, - double epsilon, - int *rx, - double *sbar, - int *violated_found); - -int GREEDY_ND_find_tight_rays(int nrows, - int nrays, - const double *f, - const double *rays, - const double *x, - const double *beta, - double epsilon, - int *tx); - -int GREEDY_ND_scale_to_ahull(int nrows, - int nrays, - const double *rays, - const int *rx, - const double *beta, - double epsilon, - const double *d, - double *alpha); - -#endif //MULTIROW_GREEDY_ND_H diff --git a/infinity/library/include/infinity/greedy-bsearch.h b/infinity/library/include/infinity/infinity-nd.h similarity index 52% rename from infinity/library/include/infinity/greedy-bsearch.h rename to infinity/library/include/infinity/infinity-nd.h index 49160cc..81753d5 100644 --- a/infinity/library/include/infinity/greedy-bsearch.h +++ b/infinity/library/include/infinity/infinity-nd.h @@ -13,21 +13,24 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef MULTIROW_GREEDY_BSEARCH_H -#define MULTIROW_GREEDY_BSEARCH_H +#ifndef MULTIROW_INFINITY_ND_H +#define MULTIROW_INFINITY_ND_H -int create_sfree_mip(int nrows, - int nrays, - const double *f, - const double *rays, - const double *bounds, - double e, - struct LP *lp); +int INFINITY_create_psi_lp(const struct ConvLFreeSet *lfree, struct LP *lp); -int GREEDY_BSEARCH_compute_bounds(int nrows, - int nrays, - const double *f, - const double *rays, - double *bounds); +int INFINITY_psi(const int nrows, + const double *q, + const double q_scale, + struct LP *lp, + double *value); -#endif //MULTIROW_GREEDY_BSEARCH_H +int INFINITY_pi(const int nrows, + const double *q, + const double q_scale, + struct LP *lp, + double *value); + +int INFINITY_ND_generate_lfree(const struct MultiRowModel *model, + struct ConvLFreeSet *lfree); + +#endif //MULTIROW_INFINITY_ND_H diff --git a/infinity/library/src/greedy-bsearch.c b/infinity/library/src/greedy-bsearch.c deleted file mode 100644 index b4bdfb2..0000000 --- a/infinity/library/src/greedy-bsearch.c +++ /dev/null @@ -1,246 +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 - -int create_sfree_mip(int nrows, - int nrays, - const double *f, - const double *rays, - const double *bounds, - double e, - struct LP *lp) -{ - int rval = 0; - - double rhs; - char sense; - - int rmatbeg = 0; - int* rmatind = 0; - double *rmatval = 0; - - rmatind = (int *) malloc((nrows + nrays) * sizeof(int)); - rmatval = (double *) malloc((nrows + nrays) * sizeof(double)); - abort_if(!rmatind, "could not allocate rmatind"); - abort_if(!rmatval, "could not allocate rmatval"); - - rval = LP_create(lp, "greedy"); - abort_if(rval, "LP_create failed"); - - // create x (basic) variables - for (int i = 0; i < nrows; i++) - { - rval = LP_new_col(lp, 0, -MILP_INFINITY, MILP_INFINITY, 'I'); - abort_if(rval, "LP_new_col failed"); - } - - // create s (non-basic) variables - for (int i = 0; i < nrays; i++) - { - rval = LP_new_col(lp, 1.0, 0, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); - } - - // add constraint \sum_{i=1}^m s_i \leq 1 - sense = 'L'; - rhs = 1.0; - - for (int i = 0; i < nrays; i++) - { - rmatind[i] = i + nrows; - rmatval[i] = 1.0; - } - - rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); - - // add constraints x_i - \sum_{j=1}^m min{e,e_j} s_j R_ji = f_i - for (int i = 0; i < nrows; i++) - { - int k = 0; - sense = 'E'; - rhs = f[i]; - - rmatind[k] = i; - rmatval[k] = 1.0; - k++; - - for (int j = 0; j < nrays; j++) - { - rmatind[k] = j + nrows; - rmatval[k] = -rays[nrows * j + i] * fmin(e, bounds[j]); - k++; - } - - rval = LP_add_rows(lp, 1, nrays + 1, &rhs, &sense, &rmatbeg, rmatind, - rmatval); - abort_if(rval, "LP_add_rows failed"); - } - -CLEANUP: - if (rmatind) free(rmatind); - if (rmatval) free(rmatval); - return rval; -} - -int GREEDY_BSEARCH_compute_bounds(int nrows, - int nrays, - const double *f, - const double *rays, - double *bounds) -{ - int rval = 0; - - struct LP lp; - double e_upper = 2 * GREEDY_BIG_E; - double e_lower = 0.0; - - int cplex_count = 0; - double cplex_time = 0; - - int iteration_count = 0; - - double *x = 0; - - x = (double *) malloc((nrays + nrows) * sizeof(double)); - abort_if(!x, "could not allocate x"); - - for (int i = 0; i < nrays; i++) - bounds[i] = GREEDY_BIG_E; - - for (int it = 0;; it++) - { - abort_if(it > 2*nrays, "stuck in an infinite loop"); - - log_verbose("Starting iteration %d...\n", it); - - iteration_count++; - - - int solution_found = 0; - int inner_count = 0; - - while (fabs(e_upper - e_lower) > GREEDY_MAX_GAP) - { - inner_count++; - - double e = (e_upper + e_lower) / 2; - log_verbose(" e=%.12lf\n", e); - - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); - - rval = create_sfree_mip(nrows, nrays, f, rays, bounds, e, &lp); - abort_if(rval, "create_sfree_mip failed"); - - if_verbose_level - { - rval = LP_write(&lp, "greedy.lp"); - abort_if(rval, "LP_write failed"); - } - - int infeasible; - cplex_count++; - - double initial_time = get_user_time(); - - log_verbose(" Optimizing...\n"); - rval = LP_optimize(&lp, &infeasible); - if (rval) - { - // Workaround for CPLEX bug. If CPLEX tell us that this problem - // is unbounded, we disable presolve and try again. - LP_free(&lp); - LP_open(&lp); - - rval = create_sfree_mip(nrows, nrays, f, rays, bounds, e, &lp); - abort_if(rval, "create_sfree_mip failed"); - - LP_disable_presolve(&lp); - - rval = LP_optimize(&lp, &infeasible); - abort_if(rval, "LP_optimize failed"); - } - - cplex_time += get_user_time() - initial_time; - - if (infeasible) - { - e_lower = e; - log_verbose(" infeasible\n"); - if (e > GREEDY_BIG_E-1) - { - LP_free(&lp); - goto OUT; - } - } - else - { - log_verbose(" feasible\n"); - e_upper = e; - solution_found = 1; - - rval = LP_get_x(&lp, x); - abort_if(rval, "LP_get_x failed"); - } - - LP_free(&lp); - } - - if (solution_found) - { - for (int j = 0; j < nrays; j++) - { - if (!DOUBLE_geq(x[nrows + j], 0.001)) continue; - bounds[j] = fmin(bounds[j] * 0.99, e_lower * 0.99); - } - } - - log_verbose(" %d iterations %12.8lf gap\n", inner_count, e_upper - - e_lower); - - e_lower = e_upper; - e_upper = 2 * GREEDY_BIG_E; - } - -OUT: - log_debug(" %6d IPs (%.2lfms per call, %.2lfs total)\n", cplex_count, - cplex_time * 1000.0 / cplex_count, cplex_time); - - for(int i = 0; i < nrays; i++) - abort_if(DOUBLE_iszero(bounds[i]), "bounds should be positive"); - - if_verbose_level - { - time_printf("Bounds:\n"); - for (int k = 0; k < nrays; k++) - time_printf(" %12.8lf %12.8lf\n", k, bounds[k], 1 / bounds[k]); - } - -CLEANUP: - if (x) free(x); - return rval; -} diff --git a/infinity/library/src/greedy-nd.c b/infinity/library/src/infinity-nd.c similarity index 86% rename from infinity/library/src/greedy-nd.c rename to infinity/library/src/infinity-nd.c index 1adfdf9..6203f38 100644 --- a/infinity/library/src/greedy-nd.c +++ b/infinity/library/src/infinity-nd.c @@ -21,8 +21,7 @@ #include #include -#include -#include +#include static long lp_count = 0; static double lp_time = 0; @@ -42,232 +41,84 @@ static double sfree_mip_time = 0; static long scale_ahull_lp_count = 0; static double scale_ahull_lp_time = 0; -static int find_interior_point_enum(const int nrows, - const int nrays, - const double *f, - const double *rays, - const double *beta, - const double epsilon, - double *x, - int *found) -{ - int rval = 0; - - int M = 1; - struct LP lp; - double best_value = INFINITY; - double *beta2 = 0; - - abort_if(nrows != 3, "not implemented"); - - beta2 = (double *) malloc(nrays * sizeof(double)); - abort_if(!beta2, "could not allocate beta2"); - - for(int i = 0; i < nrays; i++) - { - beta2[i] = fmin(epsilon, beta[i]); - } - - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); - - struct ConvLFreeSet lfree; - lfree.f = (double*) f; - lfree.nrows = nrows; - lfree.rays.dim = nrows; - lfree.rays.nrays = nrays; - lfree.rays.values = (double*) rays; - lfree.beta = beta2; - - rval = GREEDY_create_psi_lp(&lfree, &lp); - abort_if(rval, "GREEDY_create_psi_lp failed"); - - *found = 0; - - for(int x1 = -M; x1 <= M; x1++) - for(int x2 = -M; x2 <= M; x2++) - for(int x3 = -M; x3 <= M; x3++) - { - double value; - double q[3] = {x1 - f[0], - x2 - f[1], - x3 - f[2]}; - - rval = GREEDY_ND_psi(nrows, q, 1, &lp, &value); - abort_if(rval, "GREEDY_ND_psi failed"); - - if(value < best_value) - { - best_value = value; - x[0] = x1; - x[1] = x2; - x[2] = x3; - } - } - - if(best_value < 0.999) *found = 1; - -CLEANUP: - if(beta2) free(beta2); - LP_free(&lp); - return rval; -} - -static int find_interior_point_cplex(const int nrows, - const int nrays, - const double *f, - const double *rays, - const double *beta, - const double epsilon, - double *x, - int *found) +static int create_sfree_mip(int nrows, + int nrays, + const double *f, + const double *rays, + const double *bounds, + double e, + struct LP *lp) { int rval = 0; - struct LP lp; - double initial_time; - int infeasible; - double objval; - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); + double rhs; + char sense; - lp_count++; - sfree_mip_count++; - initial_time = get_user_time(); + int rmatbeg = 0; + int* rmatind = 0; + double *rmatval = 0; - rval = create_sfree_mip(nrows, nrays, f, rays, beta, epsilon, &lp); - abort_if(rval, "greate_sfree_mip failed"); + rmatind = (int *) malloc((nrows + nrays) * sizeof(int)); + rmatval = (double *) malloc((nrows + nrays) * sizeof(double)); + abort_if(!rmatind, "could not allocate rmatind"); + abort_if(!rmatval, "could not allocate rmatval"); - log_debug(" solving sfree mip...\n"); - rval = LP_optimize(&lp, &infeasible); - if(rval == ERR_MIP_TIMEOUT) goto CLEANUP; - abort_if(rval, "LP_optimize failed"); + rval = LP_create(lp, "greedy"); + abort_if(rval, "LP_create failed"); - if(infeasible) + // create x (basic) variables + for (int i = 0; i < nrows; i++) { - log_debug(" mip is infeasible. Stopping.\n"); - *found = 0; - goto CLEANUP; + rval = LP_new_col(lp, 0, -MILP_INFINITY, MILP_INFINITY, 'I'); + abort_if(rval, "LP_new_col failed"); } - rval = LP_get_x(&lp, x); - abort_if(rval, "LP_get_x failed"); - - if_verbose_level + // create s (non-basic) variables + for (int i = 0; i < nrays; i++) { - for(int i = 0; i < nrows; i++) - log_verbose(" x%d = %.8lf\n", i, x[i]); - - for(int i = 0; i < nrays; i++) - if(x[i + nrows] > 0.00001) - log_verbose(" t%d = %.8lf\n", i, x[i + nrows]); + rval = LP_new_col(lp, 1.0, 0, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col failed"); } - rval = LP_get_obj_val(&lp, &objval); - abort_if(rval, "LP_get_obj_val failed"); - - log_debug(" obj = %.8lf\n", objval); + // add constraint \sum_{i=1}^m s_i \leq 1 + sense = 'L'; + rhs = 1.0; - if(objval >= 0.999) + for (int i = 0; i < nrays; i++) { - log_debug(" set is lattice-free\n"); - *found = 0; - goto CLEANUP; + rmatind[i] = i + nrows; + rmatval[i] = 1.0; } - *found = 1; - - sfree_mip_time += get_user_time() - initial_time; - lp_time += get_user_time() - initial_time; - -CLEANUP: - LP_free(&lp); - return rval; -} - -static int bound(int nrows, - int nrays, - const double *f, - const double *rays, - const double *x, - const double *beta, - double *epsilon, - int *tx) -{ - int rval = 0; - - int found; - int *rx = 0; - double *fbar = 0; - double *sbar = 0; - - double prev_epsilon; - int count = 0; - *epsilon = GREEDY_BIG_E; - - rx = (int *) malloc(nrays * sizeof(int)); - fbar = (double *) malloc(nrows * sizeof(double)); - sbar = (double *) malloc(nrays * sizeof(double)); - - abort_if(!rx, "could not allocate rx"); - abort_if(!fbar, "could not allocate fbar"); - abort_if(!sbar, "could not allocate sbar"); + rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); - while(1) + // add constraints x_i - \sum_{j=1}^m min{e,e_j} s_j R_ji = f_i + for (int i = 0; i < nrows; i++) { - count++; - abort_if(count > 100, "infinite loop"); - - rval = GREEDY_ND_find_violated_cone(nrows, nrays, f, rays, x, beta, - *epsilon, rx, sbar, &found); - abort_if(rval, "GREEDY_ND_find_violated_cone failed"); - - if(!found) break; + int k = 0; + sense = 'E'; + rhs = f[i]; - for(int i = 0; i < nrows; i++) - fbar[i] = x[i]; + rmatind[k] = i; + rmatval[k] = 1.0; + k++; - for(int j = 0; j < nrays; j++) + for (int j = 0; j < nrays; j++) { - if(!rx[j]) continue; - const double *r = &rays[nrows * j]; - - for(int i = 0; i < nrows; i++) - fbar[i] -= min(*epsilon, beta[j]) * r[i] * sbar[j]; + rmatind[k] = j + nrows; + rmatval[k] = -rays[nrows * j + i] * fmin(e, bounds[j]); + k++; } - log_verbose("%.12lf %.12lf\n", f[0], f[1]); - log_verbose("%.12lf %.12lf\n", fbar[0], fbar[1]); - - prev_epsilon = *epsilon; - - rval = GREEDY_ND_cone_bound(nrows, nrays, fbar, rays, rx, x, beta, - epsilon); - abort_if(rval, "GREEDY_ND_cone_bound failed"); - - log_verbose(" e=%.12lf\n", *epsilon); - abort_if(prev_epsilon < *epsilon, "epsilon should never increase"); - } - - for(int i = 0; i < nrays; i++) - tx[i] = 0; - - if(DOUBLE_geq(*epsilon, GREEDY_BIG_E)) - { - *epsilon = INFINITY; - goto CLEANUP; - } - else - { - rval = GREEDY_ND_find_tight_rays(nrows, nrays, fbar, rays, x, beta, - *epsilon, tx); - abort_if(rval, "GREEDY_ND_find_tight_rays failed"); + rval = LP_add_rows(lp, 1, nrays + 1, &rhs, &sense, &rmatbeg, rmatind, + rmatval); + abort_if(rval, "LP_add_rows failed"); } CLEANUP: - if(sbar) free(sbar); - if(fbar) free(fbar); - if(rx) free(rx); + if (rmatind) free(rmatind); + if (rmatval) free(rmatval); return rval; } @@ -422,14 +273,15 @@ CLEANUP: return rval; } -static int create_violated_cone_lp(int nrows, - int nrays, - const double *f, - const double *rays, - const double *x, - const double *beta, - double epsilon, - struct LP *lp) +static int create_tight_rays_lp(int nrows, + int nrays, + const double *f, + const double *rays, + const double *x, + const double *beta, + double epsilon, + double delta, + struct LP *lp) { int rval = 0; @@ -445,17 +297,24 @@ static int create_violated_cone_lp(int nrows, abort_if(!rmatind, "could not allocate rmatind"); abort_if(!rmatval, "could not allocate rmatval"); - rval = LP_create(lp, "violated_cone"); + rval = LP_create(lp, "tight_rays"); abort_if(rval, "LP_create failed"); - // create s variables + // create lambda variables for(int i = 0; i < nrays; i++) { rval = LP_new_col(lp, 1.0, 0.0, MILP_INFINITY, 'C'); abort_if(rval, "LP_new_col failed"); } - // create constraint x = f + \sum(min{e, beta[r]} * r * s_r) + // create s variables + for(int i = 0; i < nrays; i++) + { + rval = LP_new_col(lp, 0.0, 0.0, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col failed"); + } + + // create constraint x = f + \sum_{r \in R} min{e, beta[r]} * r * s_r for(int j = 0; j < nrows; j++) { sense = 'E'; @@ -464,7 +323,7 @@ static int create_violated_cone_lp(int nrows, for(int i = 0; i < nrays; i++) { const double *ri = &rays[i * nrows]; - rmatind[i] = i; + rmatind[i] = nrays + i; rmatval[i] = min(epsilon, beta[i]) * ri[j]; if(DOUBLE_iszero(rmatval[i])) rmatval[i] = 0.0; } @@ -474,20 +333,114 @@ static int create_violated_cone_lp(int nrows, abort_if(rval, "LP_add_rows failed"); } - rval = LP_relax(lp); - abort_if(rval, "LP_relax failed"); + // create constraint \sum_{r \in R} s_r = 1 + sense = 'E'; + rhs = 1.0; - //rval = LP_write(lp, "violated-cone.lp"); - //abort_if(rval, "LP_write failed"); - //UTIL_pause(); + for(int i = 0; i < nrays; i++) + { + rmatind[i] = nrays + i; + rmatval[i] = 1.0; + } -CLEANUP: - if(rmatind) free(rmatind); - if(rmatval) free(rmatval); - return rval; -} + rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); -static int create_scale_to_ahull_lp(int nrows, + // create constraints \lambda_r + s_r \geq \delta + for(int i = 0; i < nrays; i++) + { + sense = 'G'; + rhs = delta; + + rmatind[0] = i; + rmatval[0] = 1.0; + + rmatind[1] = nrays + i; + rmatval[1] = 1.0; + + rval = LP_add_rows(lp, 1, 2, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); + } + + 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); + return rval; +} + +static int create_violated_cone_lp(int nrows, + int nrays, + const double *f, + const double *rays, + const double *x, + const double *beta, + double epsilon, + struct LP *lp) +{ + int rval = 0; + + double rhs; + char sense; + + int rmatbeg = 0; + int *rmatind = 0; + double *rmatval = 0; + + rmatind = (int *) malloc(nrays * sizeof(int)); + rmatval = (double *) malloc(nrays * sizeof(double)); + abort_if(!rmatind, "could not allocate rmatind"); + abort_if(!rmatval, "could not allocate rmatval"); + + rval = LP_create(lp, "violated_cone"); + abort_if(rval, "LP_create failed"); + + // create s variables + for(int i = 0; i < nrays; i++) + { + rval = LP_new_col(lp, 1.0, 0.0, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col failed"); + } + + // create constraint x = f + \sum(min{e, beta[r]} * r * s_r) + for(int j = 0; j < nrows; j++) + { + sense = 'E'; + rhs = x[j] - f[j]; + + for(int i = 0; i < nrays; i++) + { + const double *ri = &rays[i * nrows]; + rmatind[i] = i; + rmatval[i] = min(epsilon, beta[i]) * ri[j]; + if(DOUBLE_iszero(rmatval[i])) rmatval[i] = 0.0; + } + + rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, + rmatval); + abort_if(rval, "LP_add_rows failed"); + } + + rval = LP_relax(lp); + abort_if(rval, "LP_relax failed"); + + //rval = LP_write(lp, "violated-cone.lp"); + //abort_if(rval, "LP_write failed"); + //UTIL_pause(); + +CLEANUP: + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); + return rval; +} + +static int create_scale_to_ahull_lp(int nrows, int nrays, const double *rays, const int *rx, @@ -583,754 +536,880 @@ CLEANUP: return rval; } -static int create_tight_rays_lp(int nrows, - int nrays, - const double *f, - const double *rays, - const double *x, - const double *beta, - double epsilon, - double delta, - struct LP *lp) +static int find_interior_point_enum(const int nrows, + const int nrays, + const double *f, + const double *rays, + const double *beta, + const double epsilon, + double *x, + int *found) { int rval = 0; - double rhs; - char sense; - - int rmatbeg = 0; - int *rmatind = 0; - double *rmatval = 0; - - rmatind = (int *) malloc(nrays * sizeof(int)); - rmatval = (double *) malloc(nrays * sizeof(double)); - abort_if(!rmatind, "could not allocate rmatind"); - abort_if(!rmatval, "could not allocate rmatval"); - - rval = LP_create(lp, "tight_rays"); - abort_if(rval, "LP_create failed"); - - // create lambda variables - for(int i = 0; i < nrays; i++) - { - rval = LP_new_col(lp, 1.0, 0.0, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); - } - - // create s variables - for(int i = 0; i < nrays; i++) - { - rval = LP_new_col(lp, 0.0, 0.0, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); - } - - // create constraint x = f + \sum_{r \in R} min{e, beta[r]} * r * s_r - for(int j = 0; j < nrows; j++) - { - sense = 'E'; - rhs = x[j] - f[j]; - - for(int i = 0; i < nrays; i++) - { - const double *ri = &rays[i * nrows]; - rmatind[i] = nrays + i; - rmatval[i] = min(epsilon, beta[i]) * ri[j]; - if(DOUBLE_iszero(rmatval[i])) rmatval[i] = 0.0; - } + int M = 1; + struct LP lp; + double best_value = INFINITY; + double *beta2 = 0; - rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, - rmatval); - abort_if(rval, "LP_add_rows failed"); - } + abort_if(nrows != 3, "not implemented"); - // create constraint \sum_{r \in R} s_r = 1 - sense = 'E'; - rhs = 1.0; + beta2 = (double *) malloc(nrays * sizeof(double)); + abort_if(!beta2, "could not allocate beta2"); for(int i = 0; i < nrays; i++) { - rmatind[i] = nrays + i; - rmatval[i] = 1.0; + beta2[i] = fmin(epsilon, beta[i]); } - rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); + rval = LP_open(&lp); + abort_if(rval, "LP_open failed"); - // create constraints \lambda_r + s_r \geq \delta - for(int i = 0; i < nrays; i++) - { - sense = 'G'; - rhs = delta; + struct ConvLFreeSet lfree; + lfree.f = (double*) f; + lfree.nrows = nrows; + lfree.rays.dim = nrows; + lfree.rays.nrays = nrays; + lfree.rays.values = (double*) rays; + lfree.beta = beta2; - rmatind[0] = i; - rmatval[0] = 1.0; + rval = INFINITY_create_psi_lp(&lfree, &lp); + abort_if(rval, "INFINITY_create_psi_lp failed"); - rmatind[1] = nrays + i; - rmatval[1] = 1.0; + *found = 0; - rval = LP_add_rows(lp, 1, 2, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); - } + for(int x1 = -M; x1 <= M; x1++) + for(int x2 = -M; x2 <= M; x2++) + for(int x3 = -M; x3 <= M; x3++) + { + double value; + double q[3] = {x1 - f[0], + x2 - f[1], + x3 - f[2]}; - rval = LP_relax(lp); - abort_if(rval, "LP_relax failed"); + rval = INFINITY_psi(nrows, q, 1, &lp, &value); + abort_if(rval, "INFINITY_psi failed"); - //rval = LP_write(lp, "tight-rays.lp"); - //abort_if(rval, "LP_write failed"); + if(value < best_value) + { + best_value = value; + x[0] = x1; + x[1] = x2; + x[2] = x3; + } + } + if(best_value < 0.999) *found = 1; CLEANUP: - if(rmatind) free(rmatind); - if(rmatval) free(rmatval); + if(beta2) free(beta2); + LP_free(&lp); return rval; } -#ifndef TEST_SOURCE - -int GREEDY_create_psi_lp(const struct ConvLFreeSet *lfree, struct LP *lp) +static int find_interior_point_cplex(const int nrows, + const int nrays, + const double *f, + const double *rays, + const double *beta, + const double epsilon, + double *x, + int *found) { int rval = 0; + struct LP lp; + double initial_time; + int infeasible; + double objval; - int rmatbeg = 0; - int *rmatind = 0; - double *rmatval = 0; - - int nrows = lfree->nrows; - int nrays = lfree->rays.nrays; - double *rays = lfree->rays.values; - double *beta = lfree->beta; + rval = LP_open(&lp); + abort_if(rval, "LP_open failed"); - rmatind = (int *) malloc(nrays * sizeof(int)); - rmatval = (double *) malloc(nrays * sizeof(double)); - abort_if(!rmatind, "could not allocate rmatind"); - abort_if(!rmatval, "could not allocate rmatval"); + lp_count++; + sfree_mip_count++; + initial_time = get_user_time(); - rval = LP_create(lp, "psi"); - abort_if(rval, "LP_create failed"); + rval = create_sfree_mip(nrows, nrays, f, rays, beta, epsilon, &lp); + abort_if(rval, "greate_sfree_mip failed"); - // create lambda variables - for(int i = 0; i < nrays; i++) + 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) { - rval = LP_new_col(lp, 1.0, 0, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); + log_debug(" mip is infeasible. Stopping.\n"); + *found = 0; + goto CLEANUP; } - // create constraint 0 = \sum_{i=1}^m \lambda_i r^i_j beta_i - for(int j = 0; j < nrows; j++) + rval = LP_get_x(&lp, x); + abort_if(rval, "LP_get_x failed"); + + if_verbose_level { - char sense = 'E'; - double rhs = 0; - int nz = 0; + for(int i = 0; i < nrows; i++) + log_verbose(" x%d = %.8lf\n", i, x[i]); for(int i = 0; i < nrays; i++) - { - const double *ri = &rays[i * nrows]; - rmatind[nz] = i; - rmatval[nz] = ri[j] * beta[i]; - nz++; - } - - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); + if(x[i + nrows] > 0.00001) + log_verbose(" t%d = %.8lf\n", i, x[i + nrows]); } - rval = LP_relax(lp); - abort_if(rval, "LP_relax failed"); + rval = LP_get_obj_val(&lp, &objval); + abort_if(rval, "LP_get_obj_val failed"); - if_verbose_level + log_debug(" obj = %.8lf\n", objval); + + if(objval >= 0.999) { - rval = LP_write(lp, "psi.lp"); - abort_if(rval, "LP_write failed"); + log_debug(" set is lattice-free\n"); + *found = 0; + goto CLEANUP; } + *found = 1; + + sfree_mip_time += get_user_time() - initial_time; + lp_time += get_user_time() - initial_time; + CLEANUP: - if(rmatind) free(rmatind); - if(rmatval) free(rmatval); + LP_free(&lp); return rval; } -int GREEDY_ND_pi(const int nrows, - const double *q, - const double q_scale, - struct LP *lp, - double *value) +static int cone_bound(int nrows, + int nrays, + const double *f, + const double *rays, + const int *rx, + const double *x, + const double *beta, + double *epsilon) { + log_verbose(" find_epsilon\n"); int rval = 0; - abort_if(nrows > 3, "not implemented"); + int *t = 0; + long it = 0; - double best_value = 1e100; + t = (int *) malloc(nrays * sizeof(int)); + abort_if(!t, "could not allocate t"); - if(nrows == 2) + for(int i = 0; i < nrays; i++) + t[i] = 0; + + while(1) { - int M = 2; - for(int k0 = -M; k0 <= M; k0++) - for(int k1 = -M; k1 <= M; k1++) - { - double v; - double qk[] = {frac(q[0] * q_scale) + k0, - frac(q[1] * q_scale) + k1}; + it++; + log_verbose("Starting iteration %d...\n", it); - rval = GREEDY_ND_psi(nrows, qk, 1, lp, &v); - abort_if(rval, "GREEDY_ND_psi failed"); + struct LP lp; - best_value = min(best_value, v); - } - } + double initial_time = get_user_time(); + rval = LP_open(&lp); + abort_if(rval, "LP_open failed"); - if(nrows == 3) - { - int M = 2; - for(int k0 = -M; k0 <= M; k0++) - for(int k1 = -M; k1 <= M; k1++) - for(int k2 = -M; k2 <= M; k2++) - { - double v; - double qk[] = {frac(q[0] * q_scale) + k0, - frac(q[1] * q_scale) + k1, - frac(q[2] * q_scale) + k2}; + rval = create_find_epsilon_lp(nrows, nrays, f, rays, t, rx, x, beta, + &lp); + abort_if(rval, "create_find_epsilon_lp failed"); - rval = GREEDY_ND_psi(nrows, qk, 1, lp, &v); - abort_if(rval, "GREEDY_ND_psi failed"); + int infeasible; + rval = LP_optimize(&lp, &infeasible); + abort_if(rval, "LP_optimize failed"); - best_value = min(best_value, v); - } - } + lp_count++; + lp_time += get_user_time() - initial_time; - *value = best_value; + epsilon_lp_count++; + epsilon_lp_time += get_user_time() - initial_time; -CLEANUP: - return rval; -} + if(infeasible) + { + *epsilon = INFINITY; + log_verbose(" infeasible\n"); + LP_free(&lp); + goto CLEANUP; + } -/* - * Given a point f, a list of rays r1,...,rm, some real numbers b1,...,bm, a - * point q, and a real number q_scale, this function evaluates psi_B(q * - * q_scale), where B is the convex hull of {f + ri * bi}_i=1^m. - */ -int GREEDY_ND_psi(const int nrows, - const double *q, - const double q_scale, - struct LP *lp, - double *value) -{ - int rval = 0; + double obj; + rval = LP_get_obj_val(&lp, &obj); + abort_if(rval, "LP_get_obj_val failed"); - for(int j = 0; j < nrows; j++) - { - rval = LP_change_rhs(lp, j, q[j] * q_scale); - abort_if(rval, "LP_change_rhs failed"); - } + log_verbose(" obj=%.6lf\n", obj); - int infeasible; - rval = LP_optimize(lp, &infeasible); - abort_if(rval, "LP_optimize failed"); + for(int i = 0; i < nrays; i++) + { + if(!rx[i]) continue; + log_verbose(" beta[%d]=%.6lf\n", i, beta[i]); + } - if(infeasible) - { - *value = INFINITY; - } - else - { - rval = LP_get_obj_val(lp, value); - abort_if(rval, "LP_get_obj_val failed"); + LP_free(&lp); + + double e_min = INFINITY; + double e_max = -INFINITY; + + for(int i = 0; i < nrays; i++) + { + if(!rx[i] || t[i]) continue; + e_min = min(e_min, beta[i]); + } + + for(int i = 0; i < nrays; i++) + { + if(!rx[i]) continue; + e_max = fmax(e_max, beta[i]); + } + + log_verbose(" e_max=%.6lf\n", e_max); + log_verbose(" e_min=%.6lf\n", e_min); + + if(DOUBLE_leq(obj, e_min)) + { + if(DOUBLE_geq(obj, e_max)) + *epsilon = INFINITY; + else + *epsilon = obj; + + goto CLEANUP; + } + else + { + for(int i = 0; i < nrays; i++) + if(rx[i] && DOUBLE_eq(beta[i], e_min)) + t[i] = 1; + } } CLEANUP: + log_verbose(" e=%.6lf\n", *epsilon); + if(t) free(t); return rval; } -int INFINITY_ND_generate_lfree(const struct MultiRowModel *model, - struct ConvLFreeSet *lfree) +static int find_tight_rays(int nrows, + int nrays, + const double *f, + const double *rays, + const double *x, + const double *beta, + double epsilon, + int *tx) { + log_verbose(" find_tight_rays\n"); int rval = 0; - int nrows = model->nrows; - int nrays = model->rays.nrays; + int infeasible = 0; + double initial_time = 0; + const double delta = 0.001; - double *f = lfree->f; - double *rays = lfree->rays.values; - double *beta = lfree->beta; + struct LP lp; + double *sbar = 0; - lfree->nrows = model->nrows; - lfree->rays.nrays = nrays; - memcpy(rays, model->rays.values, nrays * nrows * sizeof(double)); - memcpy(f, model->f, nrows * sizeof(double)); + sbar = (double *) malloc(2 * nrays * sizeof(double)); + abort_if(!sbar, "could not allocate sbar"); - double *x = 0; + initial_time = get_user_time(); + rval = LP_open(&lp); + abort_if(rval, "LP_open failed"); - int *t = 0; - int *tx = 0; + rval = create_tight_rays_lp(nrows, nrays, f, rays, x, beta, epsilon, delta, + &lp); + abort_if(rval, "create_tight_rays_lp failed"); - t = (int *) malloc(nrays * sizeof(int)); - tx = (int *) malloc(nrays * sizeof(int)); - abort_if(!t, "could not allocate t"); - abort_if(!tx, "could not allocate tx"); + rval = LP_optimize(&lp, &infeasible); + abort_if(rval, "LP_optimize failed"); - x = (double *) malloc((nrows + nrays) * sizeof(double)); - abort_if(!x, "could not allocate x"); + lp_count++; + lp_time += get_user_time() - initial_time; - for(int i = 0; i < nrays; i++) - beta[i] = GREEDY_BIG_E; + tight_lp_count++; + tight_lp_time += get_user_time() - initial_time; - int it = 0; + abort_if(infeasible, "tight_rays_lp is infeasible"); - //lp_time = 0; - //lp_count = 0; + rval = LP_get_x(&lp, sbar); + abort_if(rval, "LP_get_x failed"); - //epsilon_lp_count = 0; - //epsilon_lp_time = 0; - // - //sfree_mip_count = 0; - //sfree_mip_time = 0; + for(int i = 0; i < nrays; i++) + tx[i] = DOUBLE_iszero(sbar[i]); - //tight_lp_count = 0; - //tight_lp_time = 0; + for(int i = 0; i < nrays; i++) + log_verbose(" tx[%d]=%d\n", i, tx[i]); - //violated_lp_count = 0; - //violated_lp_time = 0; +CLEANUP: + if(sbar) free(sbar); + LP_free(&lp); + return rval; +} - //scale_ahull_lp_time = 0; - //scale_ahull_lp_count = 0; +static int find_violated_cone(int nrows, + int nrays, + const double *f, + const double *rays, + const double *x, + const double *beta, + double epsilon, + int *rx, + double *sbar, + int *violated_found) +{ + log_verbose(" find_violated_cone\n"); + int rval = 0; - long x_count = 0; - double epsilon; + struct LP lp; - for(int i = 0; i < nrows; i++) - log_verbose(" f[%d] = %.12lf\n", i, f[i]); + double initial_time = get_user_time(); - while(1) - { - it++; - abort_if(it > 10 * nrays, "infinite loop"); + rval = LP_open(&lp); + abort_if(rval, "LP_open failed"); - log_debug("Starting iteration %d...\n", it); - epsilon = INFINITY; + rval = create_violated_cone_lp(nrows, nrays, f, rays, x, beta, epsilon, + &lp); + abort_if(rval, "create_violated_cone_lp failed"); - for(int i = 0; i < nrays; i++) - t[i] = 0; + int infeasible; + rval = LP_optimize(&lp, &infeasible); + abort_if(rval, "LP_optimize failed"); - for(int i = 0; i < nrows; i++) - x[i] = 0; + lp_count++; + lp_time += get_user_time() - initial_time; - while(1) - { - log_debug(" epsilon = %.8lf\n", epsilon); + violated_lp_count++; + violated_lp_time += get_user_time() - initial_time; - int found = 0; + rval = LP_get_x(&lp, sbar); + abort_if(rval, "LP_get_x failed"); - if(nrows == 3) - { - rval = find_interior_point_enum(nrows, nrays, f, - model->rays.values, beta, epsilon, x, &found); - abort_if(rval, "find_interior_point_enum failed"); - } + for(int i = 0; i < nrays; i++) + rx[i] = 0; - if(!found) - { - rval = find_interior_point_cplex(nrows, nrays, f, - model->rays.values, beta, epsilon, x, &found); - if(rval == ERR_MIP_TIMEOUT) goto CLEANUP; - abort_if(rval, "find_interior_point_cplex failed"); - if(!found) break; - } + if(infeasible) + goto CLEANUP; - log_debug(" found interior point:\n"); - for(int i = 0; i < nrows; i++) log_debug(" %.2lf\n", x[i]); + double obj; + rval = LP_get_obj_val(&lp, &obj); + abort_if(rval, "LP_get_obj_val failed"); - x_count++; - abort_if(x_count > 1000, "infinite loop"); + log_verbose(" o=%.8lf\n", obj); - double epsilon_x; - rval = bound(nrows, nrays, f, model->rays.values, x, beta, - &epsilon_x, tx); - abort_if(rval, "bound failed"); -// epsilon_x *= 0.999; + if(DOUBLE_geq(obj, 0.999)) + { + *violated_found = 0; + } + else + { + *violated_found = 1; - if(isinf(epsilon_x)) break; + log_verbose("Violated cone found\n"); + log_verbose(" f=%.8lf %.8lf\n", f[0], f[1]); + log_verbose(" x=%.8lf %.8lf\n", x[0], x[1]); - log_debug(" epsilon_x = %.8lf\n", epsilon_x); + for(int i = 0; i < nrays; i++) + { + rx[i] = (sbar[i] > 1e-9); - if(DOUBLE_eq(epsilon_x, epsilon)) + if(rx[i]) if_verbose_level { - for(int i = 0; i < nrays; i++) - if(tx[i]) t[i] = 1; - } - else if(epsilon_x < epsilon) - { - epsilon = epsilon_x; - for(int i = 0; i < nrays; i++) - t[i] = tx[i]; + double m = min(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]=%.8lf %.8lf\n", i, m * r[0], m * r[1]); } } + } - if(isinf(epsilon)) - break; +CLEANUP: + LP_free(&lp); + return rval; +} - int skip_ahull = 0; +static int bound(int nrows, + int nrays, + const double *f, + const double *rays, + const double *x, + const double *beta, + double *epsilon, + int *tx) +{ + int rval = 0; - for(int i = 0; i < nrays; i++) - { - if(t[i]) - { - beta[i] = min(beta[i], epsilon); -// beta[i] *= 0.999; - } - else if(!skip_ahull) - { - double alpha; - const double *d = LFREE_get_ray(&model->rays, i); + int found; + int *rx = 0; + double *fbar = 0; + double *sbar = 0; - rval = GREEDY_ND_scale_to_ahull(nrows, nrays, - model->rays.values, t, beta, epsilon, d, &alpha); - abort_if(rval, "GREEDY_ND_scale_to_ahull failed"); + double prev_epsilon; + int count = 0; + *epsilon = GREEDY_BIG_E; - if(DOUBLE_iszero(alpha)) - { - skip_ahull = 1; - continue; - } + rx = (int *) malloc(nrays * sizeof(int)); + fbar = (double *) malloc(nrows * sizeof(double)); + sbar = (double *) malloc(nrays * sizeof(double)); - beta[i] = min(beta[i], alpha); -// beta[i] *= 0.999; - } + abort_if(!rx, "could not allocate rx"); + abort_if(!fbar, "could not allocate fbar"); + abort_if(!sbar, "could not allocate sbar"); - log_debug(" beta[%2d] = %.4lf\n", i, beta[i]); - } + while(1) + { + count++; + abort_if(count > 100, "infinite loop"); - log_debug("epsilon = %.6lf\n", epsilon); - } + rval = find_violated_cone(nrows, nrays, f, rays, x, beta, *epsilon, rx, + sbar, &found); + abort_if(rval, "find_violated_cone failed"); - log_debug(" %6ld lattice points, %ld iterations\n", x_count, it); + if(!found) break; - log_debug(" %6ld MIPs (%.2lf ms per call, %.0lf ms total)\n", lp_count, - lp_time * 1000 / lp_count, lp_time * 1000); + for(int i = 0; i < nrows; i++) + fbar[i] = x[i]; - log_debug( - " %6ld S-free MIPs (%.2lf ms per call, %.0lf ms total)\n", - sfree_mip_count, sfree_mip_time * 1000 / sfree_mip_count, - sfree_mip_time * 1000); + for(int j = 0; j < nrays; j++) + { + if(!rx[j]) continue; + const double *r = &rays[nrows * j]; - log_debug( - " %6ld epsilon LPs (%.2lf ms per call, %.0lf ms total)\n", - epsilon_lp_count, epsilon_lp_time * 1000 / epsilon_lp_count, - epsilon_lp_time * 1000); + for(int i = 0; i < nrows; i++) + fbar[i] -= min(*epsilon, beta[j]) * r[i] * sbar[j]; + } - log_debug( - " %6ld tight-rays LPs (%.2lf ms per call, %.0lf ms total)\n", - tight_lp_count, tight_lp_time * 1000 / tight_lp_count, - tight_lp_time * 1000); + log_verbose("%.12lf %.12lf\n", f[0], f[1]); + log_verbose("%.12lf %.12lf\n", fbar[0], fbar[1]); - log_debug( - " %6ld violated-cone LPs (%.2lf ms per call, %.0lf ms total)\n", - violated_lp_count, violated_lp_time * 1000 / violated_lp_count, - violated_lp_time * 1000); + prev_epsilon = *epsilon; - log_debug( - " %6ld scale-to-ahull LPs (%.2lf ms per call, %.0lf ms total)\n", - scale_ahull_lp_count, - scale_ahull_lp_time * 1000 / scale_ahull_lp_count, - scale_ahull_lp_time * 1000); + rval = cone_bound(nrows, nrays, fbar, rays, rx, x, beta, epsilon); + abort_if(rval, "cone_bound failed"); + + log_verbose(" e=%.12lf\n", *epsilon); + abort_if(prev_epsilon < *epsilon, "epsilon should never increase"); + } + + for(int i = 0; i < nrays; i++) + tx[i] = 0; + + if(DOUBLE_geq(*epsilon, GREEDY_BIG_E)) + { + *epsilon = INFINITY; + goto CLEANUP; + } + else + { + rval = find_tight_rays(nrows, nrays, fbar, rays, x, beta, *epsilon, tx); + abort_if(rval, "find_tight_rays failed"); + } CLEANUP: - if(x) free(x); - if(t) free(t); - if(tx) free(tx); + if(sbar) free(sbar); + if(fbar) free(fbar); + if(rx) free(rx); return rval; } -int GREEDY_ND_cone_bound(int nrows, - int nrays, - const double *f, - const double *rays, - const int *rx, - const double *x, - const double *beta, - double *epsilon) +static int scale_to_ahull(int nrows, + int nrays, + const double *rays, + const int *rx, + const double *beta, + double epsilon, + const double *d, + double *alpha) { - log_verbose(" find_epsilon\n"); + log_verbose(" scale_to_ahull\n"); int rval = 0; + *alpha = INFINITY; - int *t = 0; - long it = 0; + struct LP lp; + double *x = 0; + double initial_time; - t = (int *) malloc(nrays * sizeof(int)); - abort_if(!t, "could not allocate t"); + x = (double *) malloc((nrays + 1) * sizeof(double)); + abort_if(!x, "could not allocate x"); - for(int i = 0; i < nrays; i++) - t[i] = 0; + initial_time = get_user_time(); + rval = LP_open(&lp); + abort_if(rval, "LP_open failed"); - while(1) - { - it++; - log_verbose("Starting iteration %d...\n", it); + rval = create_scale_to_ahull_lp(nrows, nrays, rays, rx, beta, epsilon, d, + &lp); + abort_if(rval, "create_scale_to_ahull_lp failed"); - struct LP lp; + int infeasible; + rval = LP_optimize(&lp, &infeasible); + abort_if(rval, "LP_optimize failed"); - double initial_time = get_user_time(); - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); + lp_count++; + lp_time += get_user_time() - initial_time; - rval = create_find_epsilon_lp(nrows, nrays, f, rays, t, rx, x, beta, - &lp); - abort_if(rval, "create_find_epsilon_lp failed"); + scale_ahull_lp_count++; + scale_ahull_lp_time += get_user_time() - initial_time; - int infeasible; - rval = LP_optimize(&lp, &infeasible); - abort_if(rval, "LP_optimize failed"); + if(infeasible) + goto CLEANUP; - lp_count++; - lp_time += get_user_time() - initial_time; + rval = LP_get_x(&lp, x); + abort_if(rval, "LP_get_x failed"); - epsilon_lp_count++; - epsilon_lp_time += get_user_time() - initial_time; + *alpha = x[0]; - if(infeasible) - { - *epsilon = INFINITY; - log_verbose(" infeasible\n"); - LP_free(&lp); - goto CLEANUP; - } +CLEANUP: + if(x) free(x); + LP_free(&lp); + return rval; +} - double obj; - rval = LP_get_obj_val(&lp, &obj); - abort_if(rval, "LP_get_obj_val failed"); +#ifndef TEST_SOURCE - log_verbose(" obj=%.6lf\n", obj); +int INFINITY_create_psi_lp(const struct ConvLFreeSet *lfree, struct LP *lp) +{ + int rval = 0; - for(int i = 0; i < nrays; i++) - { - if(!rx[i]) continue; - log_verbose(" beta[%d]=%.6lf\n", i, beta[i]); - } + int rmatbeg = 0; + int *rmatind = 0; + double *rmatval = 0; - LP_free(&lp); + int nrows = lfree->nrows; + int nrays = lfree->rays.nrays; + double *rays = lfree->rays.values; + double *beta = lfree->beta; - double e_min = INFINITY; - double e_max = -INFINITY; + rmatind = (int *) malloc(nrays * sizeof(int)); + rmatval = (double *) malloc(nrays * sizeof(double)); + abort_if(!rmatind, "could not allocate rmatind"); + abort_if(!rmatval, "could not allocate rmatval"); - for(int i = 0; i < nrays; i++) - { - if(!rx[i] || t[i]) continue; - e_min = min(e_min, beta[i]); - } + rval = LP_create(lp, "psi"); + abort_if(rval, "LP_create failed"); + + // create lambda variables + for(int i = 0; i < nrays; i++) + { + rval = LP_new_col(lp, 1.0, 0, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col failed"); + } + + // create constraint 0 = \sum_{i=1}^m \lambda_i r^i_j beta_i + for(int j = 0; j < nrows; j++) + { + char sense = 'E'; + double rhs = 0; + int nz = 0; for(int i = 0; i < nrays; i++) { - if(!rx[i]) continue; - e_max = fmax(e_max, beta[i]); + const double *ri = &rays[i * nrows]; + rmatind[nz] = i; + rmatval[nz] = ri[j] * beta[i]; + nz++; } - log_verbose(" e_max=%.6lf\n", e_max); - log_verbose(" e_min=%.6lf\n", e_min); + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); + } - if(DOUBLE_leq(obj, e_min)) - { - if(DOUBLE_geq(obj, e_max)) - *epsilon = INFINITY; - else - *epsilon = obj; + rval = LP_relax(lp); + abort_if(rval, "LP_relax failed"); - goto CLEANUP; - } - else - { - for(int i = 0; i < nrays; i++) - if(rx[i] && DOUBLE_eq(beta[i], e_min)) - t[i] = 1; - } + if_verbose_level + { + rval = LP_write(lp, "psi.lp"); + abort_if(rval, "LP_write failed"); } CLEANUP: - log_verbose(" e=%.6lf\n", *epsilon); - if(t) free(t); + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); return rval; } -int GREEDY_ND_scale_to_ahull(int nrows, - int nrays, - const double *rays, - const int *rx, - const double *beta, - double epsilon, - const double *d, - double *alpha) +int INFINITY_pi(const int nrows, + const double *q, + const double q_scale, + struct LP *lp, + double *value) { - log_verbose(" scale_to_ahull\n"); int rval = 0; - *alpha = INFINITY; - struct LP lp; - double *x = 0; + abort_if(nrows > 3, "not implemented"); - x = (double *) malloc((nrays + 1) * sizeof(double)); - abort_if(!x, "could not allocate x"); + double best_value = 1e100; - double initial_time = get_user_time(); - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); + if(nrows == 2) + { + int M = 2; + for(int k0 = -M; k0 <= M; k0++) + for(int k1 = -M; k1 <= M; k1++) + { + double v; + double qk[] = {frac(q[0] * q_scale) + k0, + frac(q[1] * q_scale) + k1}; - rval = create_scale_to_ahull_lp(nrows, nrays, rays, rx, beta, epsilon, d, - &lp); - abort_if(rval, "create_scale_to_ahull_lp failed"); + rval = INFINITY_psi(nrows, qk, 1, lp, &v); + abort_if(rval, "INFINITY_psi failed"); + + best_value = min(best_value, v); + } + } + + if(nrows == 3) + { + int M = 2; + for(int k0 = -M; k0 <= M; k0++) + for(int k1 = -M; k1 <= M; k1++) + for(int k2 = -M; k2 <= M; k2++) + { + double v; + double qk[] = {frac(q[0] * q_scale) + k0, + frac(q[1] * q_scale) + k1, + frac(q[2] * q_scale) + k2}; + + rval = INFINITY_psi(nrows, qk, 1, lp, &v); + abort_if(rval, "INFINITY_psi failed"); + + best_value = min(best_value, v); + } + } + + *value = best_value; + +CLEANUP: + return rval; +} + +/* + * Given a point f, a list of rays r1,...,rm, some real numbers b1,...,bm, a + * point q, and a real number q_scale, this function evaluates psi_B(q * + * q_scale), where B is the convex hull of {f + ri * bi}_i=1^m. + */ +int INFINITY_psi(const int nrows, + const double *q, + const double q_scale, + struct LP *lp, + double *value) +{ + int rval = 0; + + for(int j = 0; j < nrows; j++) + { + rval = LP_change_rhs(lp, j, q[j] * q_scale); + abort_if(rval, "LP_change_rhs failed"); + } int infeasible; - rval = LP_optimize(&lp, &infeasible); + rval = LP_optimize(lp, &infeasible); abort_if(rval, "LP_optimize failed"); - lp_count++; - lp_time += get_user_time() - initial_time; - - scale_ahull_lp_count++; - scale_ahull_lp_time += get_user_time() - initial_time; - if(infeasible) - goto CLEANUP; - - rval = LP_get_x(&lp, x); - abort_if(rval, "LP_get_x failed"); - - *alpha = x[0]; + { + *value = INFINITY; + } + else + { + rval = LP_get_obj_val(lp, value); + abort_if(rval, "LP_get_obj_val failed"); + } CLEANUP: - if(x) free(x); - LP_free(&lp); return rval; } -int GREEDY_ND_find_tight_rays(int nrows, - int nrays, - const double *f, - const double *rays, - const double *x, - const double *beta, - double epsilon, - int *tx) +int INFINITY_ND_generate_lfree(const struct MultiRowModel *model, + struct ConvLFreeSet *lfree) { - log_verbose(" find_tight_rays\n"); int rval = 0; - const double delta = 0.001; + int nrows = model->nrows; + int nrays = model->rays.nrays; - struct LP lp; - double *sbar = 0; + double *f = lfree->f; + double *rays = lfree->rays.values; + double *beta = lfree->beta; - sbar = (double *) malloc(2 * nrays * sizeof(double)); - abort_if(!sbar, "could not allocate sbar"); + lfree->nrows = model->nrows; + lfree->rays.nrays = nrays; + memcpy(rays, model->rays.values, nrays * nrows * sizeof(double)); + memcpy(f, model->f, nrows * sizeof(double)); - double initial_time = get_user_time(); - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); + double *x = 0; - rval = create_tight_rays_lp(nrows, nrays, f, rays, x, beta, epsilon, delta, - &lp); - abort_if(rval, "create_tight_rays_lp failed"); + int *t = 0; + int *tx = 0; - int infeasible = 0; - rval = LP_optimize(&lp, &infeasible); - abort_if(rval, "LP_optimize failed"); + t = (int *) malloc(nrays * sizeof(int)); + tx = (int *) malloc(nrays * sizeof(int)); + abort_if(!t, "could not allocate t"); + abort_if(!tx, "could not allocate tx"); - lp_count++; - lp_time += get_user_time() - initial_time; + x = (double *) malloc((nrows + nrays) * sizeof(double)); + abort_if(!x, "could not allocate x"); - tight_lp_count++; - tight_lp_time += get_user_time() - initial_time; + for(int i = 0; i < nrays; i++) + beta[i] = GREEDY_BIG_E; - abort_if(infeasible, "tight_rays_lp is infeasible"); + int it = 0; - rval = LP_get_x(&lp, sbar); - abort_if(rval, "LP_get_x failed"); + //lp_time = 0; + //lp_count = 0; - for(int i = 0; i < nrays; i++) - tx[i] = DOUBLE_iszero(sbar[i]); + //epsilon_lp_count = 0; + //epsilon_lp_time = 0; + // + //sfree_mip_count = 0; + //sfree_mip_time = 0; - for(int i = 0; i < nrays; i++) - log_verbose(" tx[%d]=%d\n", i, tx[i]); + //tight_lp_count = 0; + //tight_lp_time = 0; -CLEANUP: - if(sbar) free(sbar); - LP_free(&lp); - return rval; -} + //violated_lp_count = 0; + //violated_lp_time = 0; -int GREEDY_ND_find_violated_cone(int nrows, - int nrays, - const double *f, - const double *rays, - const double *x, - const double *beta, - double epsilon, - int *rx, - double *sbar, - int *violated_found) -{ - log_verbose(" find_violated_cone\n"); - int rval = 0; + //scale_ahull_lp_time = 0; + //scale_ahull_lp_count = 0; - struct LP lp; + long x_count = 0; + double epsilon; - double initial_time = get_user_time(); + for(int i = 0; i < nrows; i++) + log_verbose(" f[%d] = %.12lf\n", i, f[i]); - rval = LP_open(&lp); - abort_if(rval, "LP_open failed"); + while(1) + { + it++; + abort_if(it > 10 * nrays, "infinite loop"); - rval = create_violated_cone_lp(nrows, nrays, f, rays, x, beta, epsilon, - &lp); - abort_if(rval, "create_violated_cone_lp failed"); + log_debug("Starting iteration %d...\n", it); + epsilon = INFINITY; - int infeasible; - rval = LP_optimize(&lp, &infeasible); - abort_if(rval, "LP_optimize failed"); + for(int i = 0; i < nrays; i++) + t[i] = 0; - lp_count++; - lp_time += get_user_time() - initial_time; + for(int i = 0; i < nrows; i++) + x[i] = 0; - violated_lp_count++; - violated_lp_time += get_user_time() - initial_time; + while(1) + { + log_debug(" epsilon = %.8lf\n", epsilon); - rval = LP_get_x(&lp, sbar); - abort_if(rval, "LP_get_x failed"); + int found = 0; - for(int i = 0; i < nrays; i++) - rx[i] = 0; + if(nrows == 3) + { + rval = find_interior_point_enum(nrows, nrays, f, + model->rays.values, beta, epsilon, x, &found); + abort_if(rval, "find_interior_point_enum failed"); + } - if(infeasible) - goto CLEANUP; + if(!found) + { + rval = find_interior_point_cplex(nrows, nrays, f, + model->rays.values, beta, epsilon, x, &found); + if(rval == ERR_MIP_TIMEOUT) goto CLEANUP; + abort_if(rval, "find_interior_point_cplex failed"); + if(!found) break; + } - double obj; - rval = LP_get_obj_val(&lp, &obj); - abort_if(rval, "LP_get_obj_val failed"); + log_debug(" found interior point:\n"); + for(int i = 0; i < nrows; i++) log_debug(" %.2lf\n", x[i]); - log_verbose(" o=%.8lf\n", obj); + x_count++; + abort_if(x_count > 1000, "infinite loop"); - if(DOUBLE_geq(obj, 0.999)) - { - *violated_found = 0; - } - else - { - *violated_found = 1; + double epsilon_x; + rval = bound(nrows, nrays, f, model->rays.values, x, beta, + &epsilon_x, tx); + abort_if(rval, "bound failed"); +// epsilon_x *= 0.999; - log_verbose("Violated cone found\n"); + if(isinf(epsilon_x)) break; - log_verbose(" f=%.8lf %.8lf\n", f[0], f[1]); - log_verbose(" x=%.8lf %.8lf\n", x[0], x[1]); + log_debug(" epsilon_x = %.8lf\n", epsilon_x); + + if(DOUBLE_eq(epsilon_x, epsilon)) + { + for(int i = 0; i < nrays; i++) + if(tx[i]) t[i] = 1; + } + else if(epsilon_x < epsilon) + { + epsilon = epsilon_x; + for(int i = 0; i < nrays; i++) + t[i] = tx[i]; + } + } + + if(isinf(epsilon)) + break; + + int skip_ahull = 0; for(int i = 0; i < nrays; i++) { - rx[i] = (sbar[i] > 1e-9); - - if(rx[i]) + if(t[i]) { - double m = min(epsilon, beta[i]); - const double *r = &rays[nrows * i]; - log_verbose(" r[%d]=%.8lf %.8lf\n", i, r[0], r[1]); - log_verbose(" r[%d]=%.8lf %.8lf\n", i, m * r[0], m * r[1]); + beta[i] = min(beta[i], epsilon); +// beta[i] *= 0.999; + } + else if(!skip_ahull) + { + double alpha; + const double *d = LFREE_get_ray(&model->rays, i); + + rval = scale_to_ahull(nrows, nrays, model->rays.values, t, beta, + epsilon, d, &alpha); + abort_if(rval, "scale_to_ahull failed"); + + if(DOUBLE_iszero(alpha)) + { + skip_ahull = 1; + continue; + } + + beta[i] = min(beta[i], alpha); +// beta[i] *= 0.999; } + + log_debug(" beta[%2d] = %.4lf\n", i, beta[i]); } + + log_debug("epsilon = %.6lf\n", epsilon); } + log_debug(" %6ld lattice points, %ld iterations\n", x_count, it); + + log_debug(" %6ld MIPs (%.2lf ms per call, %.0lf ms total)\n", lp_count, + lp_time * 1000 / lp_count, lp_time * 1000); + + log_debug( + " %6ld S-free MIPs (%.2lf ms per call, %.0lf ms total)\n", + sfree_mip_count, sfree_mip_time * 1000 / sfree_mip_count, + sfree_mip_time * 1000); + + log_debug( + " %6ld epsilon LPs (%.2lf ms per call, %.0lf ms total)\n", + epsilon_lp_count, epsilon_lp_time * 1000 / epsilon_lp_count, + epsilon_lp_time * 1000); + + log_debug( + " %6ld tight-rays LPs (%.2lf ms per call, %.0lf ms total)\n", + tight_lp_count, tight_lp_time * 1000 / tight_lp_count, + tight_lp_time * 1000); + + log_debug( + " %6ld violated-cone LPs (%.2lf ms per call, %.0lf ms total)\n", + violated_lp_count, violated_lp_time * 1000 / violated_lp_count, + violated_lp_time * 1000); + + log_debug( + " %6ld scale-to-ahull LPs (%.2lf ms per call, %.0lf ms total)\n", + scale_ahull_lp_count, + scale_ahull_lp_time * 1000 / scale_ahull_lp_count, + scale_ahull_lp_time * 1000); + CLEANUP: - LP_free(&lp); + if(x) free(x); + if(t) free(t); + if(tx) free(tx); return rval; } diff --git a/infinity/library/src/infinity.c b/infinity/library/src/infinity.c index 6494f8c..ede040e 100644 --- a/infinity/library/src/infinity.c +++ b/infinity/library/src/infinity.c @@ -22,7 +22,7 @@ #include #include -#include +#include #include /** @@ -117,7 +117,7 @@ static int create_cut_from_lfree(const struct Tableau *tableau, rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = GREEDY_create_psi_lp(lfree, &lp); + rval = INFINITY_create_psi_lp(lfree, &lp); abort_if(rval, "create_psi_lp failed"); cut->nz = nvars; @@ -129,13 +129,13 @@ static int create_cut_from_lfree(const struct Tableau *tableau, if(ENABLE_LIFTING && type == MILP_INTEGER) { - rval = GREEDY_ND_pi(nrows, q, map->ray_scale[i], &lp, &value); - abort_if(rval, "GREEDY_ND_pi failed"); + rval = INFINITY_pi(nrows, q, map->ray_scale[i], &lp, &value); + abort_if(rval, "INFINITY_pi failed"); } else { - rval = GREEDY_ND_psi(nrows, q, map->ray_scale[i], &lp, &value); - abort_if(rval, "GREEDY_ND_psi failed"); + rval = INFINITY_psi(nrows, q, map->ray_scale[i], &lp, &value); + abort_if(rval, "INFINITY_psi failed"); } log_verbose(" psi[%4d] = %20.12lf %d\n", map->indices[i], value); diff --git a/infinity/library/tests/greedy-nd-test.cpp b/infinity/library/tests/infinity-nd-test.cpp similarity index 77% rename from infinity/library/tests/greedy-nd-test.cpp rename to infinity/library/tests/infinity-nd-test.cpp index fdd8f56..0f566ed 100644 --- a/infinity/library/tests/greedy-nd-test.cpp +++ b/infinity/library/tests/infinity-nd-test.cpp @@ -19,13 +19,12 @@ #define TEST_SOURCE extern "C" { -#include #include #include #include #include -#include -#include "../src/greedy-nd.c" +#include +#include "../src/infinity-nd.c" } int ENABLE_LIFTING = 0; @@ -34,7 +33,7 @@ int MAX_N_ROWS = 2; int SHOULD_DUMP_CUTS = 0; int DUMP_CUT_N = 0; -TEST(GreedyNDTest, find_violated_cone_test) +TEST(InfinityNDTest, find_violated_cone_test) { int rval = 0; @@ -56,9 +55,9 @@ TEST(GreedyNDTest, find_violated_cone_test) int rx[nrays]; int found; - rval = GREEDY_ND_find_violated_cone(nrows, nrays, f, rays, x, beta, - 1.0, rx, sbar, &found); - abort_if(rval, "GREEDY_ND_find_violated_cone failed"); + rval = find_violated_cone(nrows, nrays, f, rays, x, beta, 1.0, rx, sbar, + &found); + abort_if(rval, "find_violated_cone failed"); EXPECT_TRUE(found); EXPECT_FALSE(rx[0]); @@ -66,9 +65,9 @@ TEST(GreedyNDTest, find_violated_cone_test) EXPECT_TRUE(rx[2]); EXPECT_FALSE(rx[3]); - rval = GREEDY_ND_find_violated_cone(nrows, nrays, f, rays, x, beta, - 0.5, rx, sbar, &found); - abort_if(rval, "GREEDY_ND_find_violated_cone failed"); + rval = find_violated_cone(nrows, nrays, f, rays, x, beta, 0.5, rx, sbar, + &found); + abort_if(rval, "find_violated_cone failed"); EXPECT_FALSE(found); @@ -76,7 +75,7 @@ TEST(GreedyNDTest, find_violated_cone_test) if(rval) FAIL(); } -TEST(GreedyNDTest, find_tight_rays_test_1) +TEST(InfinityNDTest, find_tight_rays_test_1) { int rval = 0; @@ -97,8 +96,8 @@ TEST(GreedyNDTest, find_tight_rays_test_1) int tx[6]; - rval = GREEDY_ND_find_tight_rays(2, 6, f, rays, x, beta, epsilon, tx); - abort_if(rval, "GREEDY_ND_find_tight_rays failed"); + rval = find_tight_rays(2, 6, f, rays, x, beta, epsilon, tx); + abort_if(rval, "find_tight_rays failed"); EXPECT_TRUE(tx[0]); EXPECT_FALSE(tx[1]); @@ -111,7 +110,7 @@ TEST(GreedyNDTest, find_tight_rays_test_1) if(rval) FAIL(); } -TEST(GreedyNDTest, find_tight_rays_test_2) +TEST(InfinityNDTest, find_tight_rays_test_2) { int rval = 0; @@ -131,8 +130,8 @@ TEST(GreedyNDTest, find_tight_rays_test_2) int tx[6]; - rval = GREEDY_ND_find_tight_rays(2, 6, f, rays, x, beta, epsilon, tx); - abort_if(rval, "GREEDY_ND_find_tight_rays failed"); + rval = find_tight_rays(2, 6, f, rays, x, beta, epsilon, tx); + abort_if(rval, "find_tight_rays failed"); EXPECT_TRUE(tx[0]); EXPECT_FALSE(tx[1]); @@ -145,7 +144,7 @@ TEST(GreedyNDTest, find_tight_rays_test_2) if(rval) FAIL(); } -TEST(GreedyNDTest, cone_bound_test_1) +TEST(InfinityNDTest, cone_bound_test_1) { int rval = 0; @@ -168,19 +167,19 @@ TEST(GreedyNDTest, cone_bound_test_1) double epsilon; - rval = GREEDY_ND_cone_bound(2, 6, f, rays, rx1, x, beta, &epsilon); - abort_if(rval, "GREEDY_ND_cone_bound failed"); + rval = cone_bound(2, 6, f, rays, rx1, x, beta, &epsilon); + abort_if(rval, "cone_bound failed"); EXPECT_NEAR(0.5, epsilon, 1e-6); - rval = GREEDY_ND_cone_bound(2, 6, f, rays, rx2, x, beta, &epsilon); - abort_if(rval, "GREEDY_ND_cone_bound failed"); + rval = cone_bound(2, 6, f, rays, rx2, x, beta, &epsilon); + abort_if(rval, "cone_bound failed"); EXPECT_NEAR(1.0, epsilon, 1e-6); CLEANUP: if(rval) FAIL(); } -TEST(GreedyNDTest, cone_bound_test_2) +TEST(InfinityNDTest, cone_bound_test_2) { int rval = 0; @@ -201,27 +200,27 @@ TEST(GreedyNDTest, cone_bound_test_2) double epsilon; - rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x1, beta1, &epsilon); - abort_if(rval, "GREEDY_ND_cone_bound failed"); + rval = cone_bound(2, 2, f, rays, rx, x1, beta1, &epsilon); + abort_if(rval, "cone_bound failed"); EXPECT_NEAR(1.0, epsilon, 1e-6); - rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x1, beta2, &epsilon); - abort_if(rval, "GREEDY_ND_cone_bound failed"); + rval = cone_bound(2, 2, f, rays, rx, x1, beta2, &epsilon); + abort_if(rval, "cone_bound failed"); EXPECT_EQ(INFINITY, epsilon); - rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x2, beta2, &epsilon); - abort_if(rval, "GREEDY_ND_cone_bound failed"); + rval = cone_bound(2, 2, f, rays, rx, x2, beta2, &epsilon); + abort_if(rval, "cone_bound failed"); EXPECT_NEAR(1.0, epsilon, 1e-6); - rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x2, beta3, &epsilon); - abort_if(rval, "GREEDY_ND_cone_bound failed"); + rval = cone_bound(2, 2, f, rays, rx, x2, beta3, &epsilon); + abort_if(rval, "cone_bound failed"); EXPECT_EQ(INFINITY, epsilon); CLEANUP: if(rval) FAIL(); } -TEST(GreedyNDTest, bound_test_1) +TEST(InfinityNDTest, bound_test_1) { int rval = 0; @@ -277,7 +276,7 @@ CLEANUP: if(rval) FAIL(); } -TEST(GreedyNDTest, psi_test) +TEST(InfinityNDTest, psi_test) { int rval = 0; @@ -309,14 +308,14 @@ TEST(GreedyNDTest, psi_test) rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = GREEDY_create_psi_lp(&lfree, &lp); - abort_if(rval, "GREEDY_create_psi_lp failed"); + rval = INFINITY_create_psi_lp(&lfree, &lp); + abort_if(rval, "INFINITY_create_psi_lp failed"); - rval = GREEDY_ND_psi(2, q1, 1.0, &lp, &value); + rval = INFINITY_psi(2, q1, 1.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); EXPECT_NEAR(value, 2.0, 1e-6); - rval = GREEDY_ND_psi(2, q2, 2.0, &lp, &value); + rval = INFINITY_psi(2, q2, 2.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); EXPECT_NEAR(value, 8.0, 1e-6); @@ -325,7 +324,7 @@ CLEANUP: if(rval) FAIL(); } -TEST(GreedyNDTest, psi_test_2) +TEST(InfinityNDTest, psi_test_2) { int rval = 0; @@ -361,14 +360,14 @@ TEST(GreedyNDTest, psi_test_2) rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = GREEDY_create_psi_lp(&lfree, &lp); - abort_if(rval, "GREEDY_create_psi_lp failed"); + rval = INFINITY_create_psi_lp(&lfree, &lp); + abort_if(rval, "INFINITY_create_psi_lp failed"); - rval = GREEDY_ND_psi(3, q1, 1.0, &lp, &value); + rval = INFINITY_psi(3, q1, 1.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); EXPECT_NEAR(value, 1.0, 1e-6); - rval = GREEDY_ND_psi(3, q2, 1.0, &lp, &value); + rval = INFINITY_psi(3, q2, 1.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); EXPECT_NEAR(value, 2.0, 1e-6); @@ -376,7 +375,7 @@ CLEANUP: if(rval) FAIL(); } -TEST(GreedyNDTest, generate_cut_test_1) +TEST(InfinityNDTest, generate_cut_test_1) { int rval = 0; @@ -415,7 +414,7 @@ CLEANUP: if(rval) FAIL(); } -TEST(GreedyNDTest, generate_cut_test_2) +TEST(InfinityNDTest, generate_cut_test_2) { int rval = 0; @@ -456,7 +455,7 @@ CLEANUP: if(rval) FAIL(); } -TEST(GreedyNTest, scale_to_ahull_test) +TEST(InfinityNDTest, scale_to_ahull_test) { int rval = 0; @@ -478,16 +477,16 @@ TEST(GreedyNTest, scale_to_ahull_test) double alpha; - rval = GREEDY_ND_scale_to_ahull(3, 4, rays, rx, beta, epsilon, d1, &alpha); - abort_if(rval, "GREEDY_ND_scale_to_ahull failed"); + rval = scale_to_ahull(3, 4, rays, rx, beta, epsilon, d1, &alpha); + abort_if(rval, "scale_to_ahull failed"); EXPECT_DOUBLE_EQ(1 / 3.0, alpha); - rval = GREEDY_ND_scale_to_ahull(3, 4, rays, rx, beta, epsilon, d2, &alpha); - abort_if(rval, "GREEDY_ND_scale_to_ahull failed"); + rval = scale_to_ahull(3, 4, rays, rx, beta, epsilon, d2, &alpha); + abort_if(rval, "scale_to_ahull failed"); EXPECT_DOUBLE_EQ(0.25, alpha); - rval = GREEDY_ND_scale_to_ahull(3, 4, rays, rx, beta, epsilon, d3, &alpha); - abort_if(rval, "GREEDY_ND_scale_to_ahull failed"); + rval = scale_to_ahull(3, 4, rays, rx, beta, epsilon, d3, &alpha); + abort_if(rval, "scale_to_ahull failed"); EXPECT_EQ(INFINITY, alpha); CLEANUP: