From 1e45817d76245baa19a4b7c40733966df4200423 Mon Sep 17 00:00:00 2001 From: Alinson Xavier Date: Sun, 30 Apr 2017 09:12:25 -0400 Subject: [PATCH] Create struct ConvLFreeSet --- infinity/library/include/infinity/greedy-nd.h | 18 +- .../library/include/infinity/infinity-2d.h | 3 +- infinity/library/src/greedy-nd.c | 1562 +++++++++-------- infinity/library/src/infinity-2d.c | 51 +- infinity/library/src/infinity.c | 82 +- infinity/library/tests/greedy-nd-test.cpp | 57 +- infinity/library/tests/infinity-2d-test.cpp | 56 +- multirow/include/multirow/cg.h | 2 +- multirow/include/multirow/lfree2d.h | 13 + multirow/include/multirow/params.h | 4 +- multirow/src/cg.c | 55 +- multirow/src/lfree2d.c | 25 + multirow/src/lp.c | 2 - multirow/tests/cg-test.cpp | 3 +- 14 files changed, 988 insertions(+), 945 deletions(-) diff --git a/infinity/library/include/infinity/greedy-nd.h b/infinity/library/include/infinity/greedy-nd.h index f948f14..0c3ec8c 100644 --- a/infinity/library/include/infinity/greedy-nd.h +++ b/infinity/library/include/infinity/greedy-nd.h @@ -16,12 +16,6 @@ #ifndef MULTIROW_GREEDY_ND_H #define MULTIROW_GREEDY_ND_H -int GREEDY_ND_next_lattice_point(int dim, - const double *lb, - const double *ub, - double *p, - int *finished); - int GREEDY_create_psi_lp(const int nrows, const int nrays, const double *f, @@ -49,16 +43,8 @@ int GREEDY_ND_pi(const int nrows, struct LP *lp, double *value); -int GREEDY_ND_generate_cut(const struct MultiRowModel *model, double *beta); - -int GREEDY_ND_bound(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); int GREEDY_ND_cone_bound(int nrows, int nrays, diff --git a/infinity/library/include/infinity/infinity-2d.h b/infinity/library/include/infinity/infinity-2d.h index ee4bc96..d0bdf75 100644 --- a/infinity/library/include/infinity/infinity-2d.h +++ b/infinity/library/include/infinity/infinity-2d.h @@ -19,6 +19,7 @@ #include -int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds); +int INFINITY_2D_generate_lfree(const struct MultiRowModel *model, + struct ConvLFreeSet *lfree); #endif //MULTIROW_INFINITY_2D_H diff --git a/infinity/library/src/greedy-nd.c b/infinity/library/src/greedy-nd.c index bbca1d9..e21690f 100644 --- a/infinity/library/src/greedy-nd.c +++ b/infinity/library/src/greedy-nd.c @@ -22,6 +22,7 @@ #include #include +#include static long lp_count = 0; static double lp_time = 0; @@ -41,14 +42,14 @@ static double sfree_mip_time = 0; static long scale_ahull_lp_count = 0; static double scale_ahull_lp_time = 0; -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) +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; @@ -105,14 +106,14 @@ CLEANUP: return rval; } -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 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; @@ -177,180 +178,692 @@ CLEANUP: return rval; } -int GREEDY_create_psi_lp(const int nrows, - const int nrays, - const double *f, - const double *rays, - const double *beta, - struct LP *lp) +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 rmatbeg = 0; - int *rmatind = 0; - double *rmatval = 0; + int found; + int *rx = 0; + double *fbar = 0; + double *sbar = 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"); + rx = (int *) malloc(nrays * sizeof(int)); + fbar = (double *) malloc(nrows * sizeof(double)); + sbar = (double *) malloc(nrays * sizeof(double)); - rval = LP_create(lp, "psi"); - abort_if(rval, "LP_create failed"); + abort_if(!rx, "could not allocate rx"); + abort_if(!fbar, "could not allocate fbar"); + abort_if(!sbar, "could not allocate sbar"); - // 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"); - } + *epsilon = GREEDY_BIG_E; + double prev_epsilon; + int count = 0; - // create constraint 0 = \sum_{i=1}^m \lambda_i r^i_j beta_i - for(int j = 0; j < nrows; j++) + while(1) { - char sense = 'E'; - double rhs = 0; - int nz = 0; + count++; + abort_if(count > 100, "infinite loop"); - for(int i = 0; i < nrays; i++) + 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; + + for(int i = 0; i < nrows; i++) + fbar[i] = x[i]; + + for(int j = 0; j < nrays; j++) { - const double *ri = &rays[i * nrows]; - rmatind[nz] = i; - rmatval[nz] = ri[j] * beta[i]; - nz++; + 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]; } - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); + 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"); } - rval = LP_relax(lp); - abort_if(rval, "LP_relax failed"); + for(int i = 0; i < nrays; i++) + tx[i] = 0; - if_verbose_level + if(DOUBLE_geq(*epsilon, GREEDY_BIG_E)) { - rval = LP_write(lp, "psi.lp"); - abort_if(rval, "LP_write failed"); + *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"); } CLEANUP: - if(rmatind) free(rmatind); - if(rmatval) free(rmatval); + if(sbar) free(sbar); + if(fbar) free(fbar); + if(rx) free(rx); return rval; } -int GREEDY_ND_pi(const int nrows, - const int nrays, - const double *f, - const double *rays, - const double *beta, - const double *q, - const double q_scale, - struct LP *lp, - double *value) +static int create_find_epsilon_lp(int nrows, + int nrays, + const double *f, + const double *rays, + const int *t, + const int *rx, + const double *x, + const double *beta, + struct LP *lp) { int rval = 0; - abort_if(nrows > 3, "not implemented"); - - double best_value = 1e100; + double rhs; + char sense; - if(nrows == 2) - { - int M = 2; - for(int k0 = -M; k0 <= M; k0++) - for(int k1 = -M; k1 <= M; k1++) - { - double qk[] = {frac(q[0] * q_scale) + k0, - frac(q[1] * q_scale) + k1}; - double value; + int nz = 0; + int *map = 0; + int rmatbeg = 0; + int *rmatind = 0; + double *rmatval = 0; - rval = GREEDY_ND_psi(nrows, nrays, f, rays, beta, qk, 1, lp, - &value); - abort_if(rval, "GREEDY_ND_psi failed"); + map = (int *) malloc(nrays * sizeof(int)); + abort_if(!map, "could not allocate map"); - best_value = min(best_value, value); - } - } + rmatind = (int *) malloc((nrays + 1 + nrows) * sizeof(int)); + rmatval = (double *) malloc((nrays + 1 + nrows) * sizeof(double)); + abort_if(!rmatind, "could not allocate rmatind"); + abort_if(!rmatval, "could not allocate rmatval"); - 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 qk[] = {frac(q[0] * q_scale) + k0, - frac(q[1] * q_scale) + k1, - frac(q[2] * q_scale) + k2}; - double value; + rval = LP_create(lp, "find_epsilon"); + abort_if(rval, "LP_create failed"); - rval = GREEDY_ND_psi(nrows, nrays, f, rays, beta, qk, 1, lp, - &value); - abort_if(rval, "GREEDY_ND_psi failed"); + // create lambda variables + int rx_count = 0; - best_value = min(best_value, value); - } - } + for(int i = 0; i < nrays + 1; i++) + { + if(i < nrays && !rx[i]) continue; - *value = best_value; + double pi = 0.0; + double lb = -MILP_INFINITY; -CLEANUP: - return rval; -} + if(i < nrays && !t[i]) + { + pi = 1.0; + lb = 0.0; + } -/* - * 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 int nrays, - const double *f, - const double *rays, - const double *beta, - const double *q, - const double q_scale, - struct LP *lp, - double *value) -{ - int rval = 0; + rval = LP_new_col(lp, pi, lb, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col 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"); + if(i < nrays) + map[i] = rx_count++; } - int infeasible; - rval = LP_optimize(lp, &infeasible); - abort_if(rval, "LP_optimize failed"); + log_verbose("rx_count=%d\n", rx_count); - if(infeasible) - { - *value = INFINITY; - } - else + // create y variables + for(int i = 0; i < nrows; i++) { - rval = LP_get_obj_val(lp, value); - abort_if(rval, "LP_get_obj_val failed"); + rval = LP_new_col(lp, 0.0, -MILP_INFINITY, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col failed"); } -CLEANUP: - return rval; -} + // create constraint y = \lambda_x x + \sum_{t \in T} \lambda_r (f + \beta_r r) + for(int j = 0; j < nrows; j++) + { + sense = 'E'; + rhs = 0.0; + nz = 0; -int GREEDY_ND_generate_cut(const struct MultiRowModel *model, double *beta) -{ - int rval = 0; - int nrows = model->nrows; - int nrays = model->rays.nrays; - double *f = model->f; + for(int i = 0; i < nrays; i++) + { + if(!t[i]) continue; + const double *ri = &rays[i * nrows]; - double *x = 0; + rmatind[nz] = map[i]; + rmatval[nz] = f[j] + beta[i] * ri[j]; + nz++; + } - int *t = 0; - int *tx = 0; + rmatind[nz] = rx_count; + rmatval[nz] = x[j]; + nz++; + + rmatind[nz] = rx_count + j + 1; + rmatval[nz] = -1.0; + nz++; + + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); + } + + // create constraint y = f + \sum_{r \in Rx \setminus T) \lambda_r r + for(int j = 0; j < nrows; j++) + { + sense = 'E'; + rhs = f[j]; + nz = 0; + + for(int i = 0; i < nrays; i++) + { + if(!rx[i] || t[i]) continue; + const double *ri = &rays[i * nrows]; + + rmatind[nz] = map[i]; + rmatval[nz] = -ri[j]; + nz++; + } + + rmatind[nz] = rx_count + j + 1; + rmatval[nz] = 1.0; + nz++; + + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); + } + + // create constraint \sum_{r \in T} \lambda_r + \lambda_x = 1 + sense = 'E'; + rhs = 1.0; + nz = 0; + + for(int i = 0; i < nrays; i++) + { + if(!t[i]) continue; + rmatind[nz] = map[i]; + rmatval[nz] = 1.0; + nz++; + } + + rmatind[nz] = rx_count; + rmatval[nz] = 1.0; + nz++; + + rval = LP_add_rows(lp, 1, nz, &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, "find-epsilon.lp"); + //abort_if(rval, "LP_write failed"); + + +CLEANUP: + if(map) free(map); + 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, + const double *beta, + double epsilon, + const double *d, + struct LP *lp) +{ + int rval = 0; + + double rhs; + char sense; + int nz; + + int rmatbeg = 0; + int *rmatind = 0; + double *rmatval = 0; + + rmatind = (int *) malloc((nrays + 1) * sizeof(int)); + rmatval = (double *) malloc((nrays + 1) * sizeof(double)); + abort_if(!rmatind, "could not allocate rmatind"); + abort_if(!rmatval, "could not allocate rmatval"); + + rval = LP_create(lp, "scale_to_ahull"); + abort_if(rval, "LP_create failed"); + + // create alpha variable + rval = LP_new_col(lp, 1.0, 0.0, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col failed"); + + + // create lambda variables + for(int i = 0; i < nrays; i++) + { + rval = LP_new_col(lp, 0.0, -MILP_INFINITY, MILP_INFINITY, 'C'); + abort_if(rval, "LP_new_col failed"); + } + + // create constraint \sum_{r \in R_x} min(e, beta[r]) * r * \lambda_r = \alpha * d + for(int j = 0; j < nrows; j++) + { + sense = 'E'; + rhs = 0.0; + nz = 0; + + rmatind[nz] = 0; + rmatval[nz] = d[j]; + nz++; + + for(int i = 0; i < nrays; i++) + { + if(!rx[i]) continue; + + const double *ri = &rays[i * nrows]; + rmatind[nz] = 1 + i; + rmatval[nz] = -min(epsilon, beta[i]) * ri[j]; + if(DOUBLE_iszero(rmatval[nz])) continue; + + nz++; + } + + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); + } + + // create constraint \sum_{r \in R_x} \lambda_r = 1 + sense = 'E'; + rhs = 1.0; + nz = 0; + + for(int i = 0; i < nrays; i++) + { + if(!rx[i]) continue; + + rmatind[nz] = 1 + i; + rmatval[nz] = 1.0; + nz++; + } + + rval = LP_add_rows(lp, 1, nz, &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, "scale-to-ahull.lp"); + //abort_if(rval, "LP_write failed"); + //UTIL_pause(); + +CLEANUP: + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); + 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) +{ + 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; + } + + rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, + rmatval); + abort_if(rval, "LP_add_rows failed"); + } + + // create constraint \sum_{r \in R} s_r = 1 + sense = 'E'; + rhs = 1.0; + + for(int i = 0; i < nrays; i++) + { + rmatind[i] = nrays + i; + rmatval[i] = 1.0; + } + + rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); + + // 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; +} + +#ifndef TEST_SOURCE + +int GREEDY_create_psi_lp(const int nrows, + const int nrays, + const double *f, + const double *rays, + const double *beta, + struct LP *lp) +{ + int rval = 0; + + 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, "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++) + { + 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"); + } + + rval = LP_relax(lp); + abort_if(rval, "LP_relax failed"); + + if_verbose_level + { + rval = LP_write(lp, "psi.lp"); + abort_if(rval, "LP_write failed"); + } + +CLEANUP: + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); + return rval; +} + +int GREEDY_ND_pi(const int nrows, + const int nrays, + const double *f, + const double *rays, + const double *beta, + const double *q, + const double q_scale, + struct LP *lp, + double *value) +{ + int rval = 0; + + abort_if(nrows > 3, "not implemented"); + + double best_value = 1e100; + + if(nrows == 2) + { + int M = 2; + for(int k0 = -M; k0 <= M; k0++) + for(int k1 = -M; k1 <= M; k1++) + { + double qk[] = {frac(q[0] * q_scale) + k0, + frac(q[1] * q_scale) + k1}; + double value; + + rval = GREEDY_ND_psi(nrows, nrays, f, rays, beta, qk, 1, lp, + &value); + abort_if(rval, "GREEDY_ND_psi failed"); + + best_value = min(best_value, value); + } + } + + 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 qk[] = {frac(q[0] * q_scale) + k0, + frac(q[1] * q_scale) + k1, + frac(q[2] * q_scale) + k2}; + double value; + + rval = GREEDY_ND_psi(nrows, nrays, f, rays, beta, qk, 1, lp, + &value); + abort_if(rval, "GREEDY_ND_psi failed"); + + best_value = min(best_value, value); + } + } + + *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 GREEDY_ND_psi(const int nrows, + const int nrays, + const double *f, + const double *rays, + const double *beta, + 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); + abort_if(rval, "LP_optimize failed"); + + if(infeasible) + { + *value = INFINITY; + } + else + { + rval = LP_get_obj_val(lp, value); + abort_if(rval, "LP_get_obj_val failed"); + } + +CLEANUP: + return rval; +} + +int INFINITY_ND_generate_lfree(const struct MultiRowModel *model, + struct ConvLFreeSet *lfree) +{ + int rval = 0; + int nrows = model->nrows; + int nrays = model->rays.nrays; + + double *f = lfree->f; + double *rays = lfree->rays.values; + double *beta = lfree->beta; + + 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 *x = 0; + + int *t = 0; + int *tx = 0; t = (int *) malloc(nrays * sizeof(int)); tx = (int *) malloc(nrays * sizeof(int)); @@ -423,347 +936,109 @@ int GREEDY_ND_generate_cut(const struct MultiRowModel *model, double *beta) if(rval == ERR_MIP_TIMEOUT) goto CLEANUP; abort_if(rval, "find_interior_point_cplex failed"); if(!found) break; - } - - log_debug(" found interior point:\n"); - for(int i = 0; i < nrows; i++) log_debug(" %.2lf\n", x[i]); - - x_count++; - abort_if(x_count > 1000, "infinite loop"); - - double epsilon_x; - rval = GREEDY_ND_bound(nrows, nrays, f, model->rays.values, x, beta, - &epsilon_x, tx); - abort_if(rval, "GREEDY_ND_bound failed"); -// epsilon_x *= 0.999; - - if(isinf(epsilon_x)) break; - - 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++) - { - 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); - - 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"); - - 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: - if(x) free(x); - if(t) free(t); - if(tx) free(tx); - return rval; -} - -int GREEDY_ND_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; - - 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"); - - *epsilon = GREEDY_BIG_E; - double prev_epsilon; - int count = 0; - - while(1) - { - 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; - - for(int i = 0; i < nrows; i++) - fbar[i] = x[i]; - - 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]; - } - - 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"); - } - -CLEANUP: - if(sbar) free(sbar); - if(fbar) free(fbar); - if(rx) free(rx); - return rval; -} - -static int create_find_epsilon_lp(int nrows, - int nrays, - const double *f, - const double *rays, - const int *t, - const int *rx, - const double *x, - const double *beta, - struct LP *lp) -{ - int rval = 0; - - double rhs; - char sense; - - int nz = 0; - int *map = 0; - int rmatbeg = 0; - int *rmatind = 0; - double *rmatval = 0; - - map = (int *) malloc(nrays * sizeof(int)); - abort_if(!map, "could not allocate map"); - - rmatind = (int *) malloc((nrays + 1 + nrows) * sizeof(int)); - rmatval = (double *) malloc((nrays + 1 + nrows) * sizeof(double)); - abort_if(!rmatind, "could not allocate rmatind"); - abort_if(!rmatval, "could not allocate rmatval"); - - rval = LP_create(lp, "find_epsilon"); - abort_if(rval, "LP_create failed"); - - // create lambda variables - int rx_count = 0; - - for(int i = 0; i < nrays + 1; i++) - { - if(i < nrays && !rx[i]) continue; - - double pi = 0.0; - double lb = -MILP_INFINITY; - - if(i < nrays && !t[i]) - { - pi = 1.0; - lb = 0.0; - } - - rval = LP_new_col(lp, pi, lb, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); + } - if(i < nrays) - map[i] = rx_count++; - } + log_debug(" found interior point:\n"); + for(int i = 0; i < nrows; i++) log_debug(" %.2lf\n", x[i]); - log_verbose("rx_count=%d\n", rx_count); + x_count++; + abort_if(x_count > 1000, "infinite loop"); - // create y variables - for(int i = 0; i < nrows; i++) - { - rval = LP_new_col(lp, 0.0, -MILP_INFINITY, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); - } + 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; - // create constraint y = \lambda_x x + \sum_{t \in T} \lambda_r (f + \beta_r r) - for(int j = 0; j < nrows; j++) - { - sense = 'E'; - rhs = 0.0; - nz = 0; + if(isinf(epsilon_x)) break; - for(int i = 0; i < nrays; i++) - { - if(!t[i]) continue; - const double *ri = &rays[i * nrows]; + log_debug(" epsilon_x = %.8lf\n", epsilon_x); - rmatind[nz] = map[i]; - rmatval[nz] = f[j] + beta[i] * ri[j]; - nz++; + 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]; + } } - rmatind[nz] = rx_count; - rmatval[nz] = x[j]; - nz++; - - rmatind[nz] = rx_count + j + 1; - rmatval[nz] = -1.0; - nz++; - - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); - } + if(isinf(epsilon)) + break; - // create constraint y = f + \sum_{r \in Rx \setminus T) \lambda_r r - for(int j = 0; j < nrows; j++) - { - sense = 'E'; - rhs = f[j]; - nz = 0; + int skip_ahull = 0; for(int i = 0; i < nrays; i++) { - if(!rx[i] || t[i]) continue; - const double *ri = &rays[i * nrows]; + 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); - rmatind[nz] = map[i]; - rmatval[nz] = -ri[j]; - nz++; - } + 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"); - rmatind[nz] = rx_count + j + 1; - rmatval[nz] = 1.0; - nz++; + if(DOUBLE_iszero(alpha)) + { + skip_ahull = 1; + continue; + } - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); - } + beta[i] = min(beta[i], alpha); +// beta[i] *= 0.999; + } - // create constraint \sum_{r \in T} \lambda_r + \lambda_x = 1 - sense = 'E'; - rhs = 1.0; - nz = 0; + log_debug(" beta[%2d] = %.4lf\n", i, beta[i]); + } - for(int i = 0; i < nrays; i++) - { - if(!t[i]) continue; - rmatind[nz] = map[i]; - rmatval[nz] = 1.0; - nz++; + log_debug("epsilon = %.6lf\n", epsilon); } - rmatind[nz] = rx_count; - rmatval[nz] = 1.0; - nz++; + log_debug(" %6ld lattice points, %ld iterations\n", x_count, it); - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); + log_debug(" %6ld MIPs (%.2lf ms per call, %.0lf ms total)\n", lp_count, + lp_time * 1000 / lp_count, lp_time * 1000); - rval = LP_relax(lp); - abort_if(rval, "LP_relax failed"); + 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); - //rval = LP_write(lp, "find-epsilon.lp"); - //abort_if(rval, "LP_write failed"); + 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: - if(map) free(map); - if(rmatind) free(rmatind); - if(rmatval) free(rmatval); + if(x) free(x); + if(t) free(t); + if(tx) free(tx); return rval; } @@ -808,167 +1083,71 @@ int GREEDY_ND_cone_bound(int nrows, abort_if(rval, "LP_optimize failed"); lp_count++; - lp_time += get_user_time() - initial_time; - - epsilon_lp_count++; - epsilon_lp_time += get_user_time() - initial_time; - - if(infeasible) - { - *epsilon = INFINITY; - log_verbose(" infeasible\n"); - LP_free(&lp); - goto CLEANUP; - } - - double obj; - rval = LP_get_obj_val(&lp, &obj); - abort_if(rval, "LP_get_obj_val failed"); - - log_verbose(" obj=%.6lf\n", obj); - - for(int i = 0; i < nrays; i++) - { - if(!rx[i]) continue; - log_verbose(" beta[%d]=%.6lf\n", i, beta[i]); - } - - 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; -} - -static int create_scale_to_ahull_lp(int nrows, - int nrays, - const double *rays, - const int *rx, - const double *beta, - double epsilon, - const double *d, - struct LP *lp) -{ - int rval = 0; - - double rhs; - char sense; - int nz; - - int rmatbeg = 0; - int *rmatind = 0; - double *rmatval = 0; - - rmatind = (int *) malloc((nrays + 1) * sizeof(int)); - rmatval = (double *) malloc((nrays + 1) * sizeof(double)); - abort_if(!rmatind, "could not allocate rmatind"); - abort_if(!rmatval, "could not allocate rmatval"); - - rval = LP_create(lp, "scale_to_ahull"); - abort_if(rval, "LP_create failed"); - - // create alpha variable - rval = LP_new_col(lp, 1.0, 0.0, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); - - - // create lambda variables - for(int i = 0; i < nrays; i++) - { - rval = LP_new_col(lp, 0.0, -MILP_INFINITY, MILP_INFINITY, 'C'); - abort_if(rval, "LP_new_col failed"); - } + lp_time += get_user_time() - initial_time; - // create constraint \sum_{r \in R_x} min(e, beta[r]) * r * \lambda_r = \alpha * d - for(int j = 0; j < nrows; j++) - { - sense = 'E'; - rhs = 0.0; - nz = 0; + epsilon_lp_count++; + epsilon_lp_time += get_user_time() - initial_time; - rmatind[nz] = 0; - rmatval[nz] = d[j]; - nz++; + if(infeasible) + { + *epsilon = INFINITY; + log_verbose(" infeasible\n"); + LP_free(&lp); + goto CLEANUP; + } + + double obj; + rval = LP_get_obj_val(&lp, &obj); + abort_if(rval, "LP_get_obj_val failed"); + + log_verbose(" obj=%.6lf\n", obj); for(int i = 0; i < nrays; i++) { if(!rx[i]) continue; - - const double *ri = &rays[i * nrows]; - rmatind[nz] = 1 + i; - rmatval[nz] = -min(epsilon, beta[i]) * ri[j]; - if(DOUBLE_iszero(rmatval[nz])) continue; - - nz++; + log_verbose(" beta[%d]=%.6lf\n", i, beta[i]); } - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); - } + LP_free(&lp); - // create constraint \sum_{r \in R_x} \lambda_r = 1 - sense = 'E'; - rhs = 1.0; - nz = 0; + double e_min = INFINITY; + double e_max = -INFINITY; - for(int i = 0; i < nrays; i++) - { - if(!rx[i]) continue; + for(int i = 0; i < nrays; i++) + { + if(!rx[i] || t[i]) continue; + e_min = min(e_min, beta[i]); + } - rmatind[nz] = 1 + i; - rmatval[nz] = 1.0; - nz++; - } + for(int i = 0; i < nrays; i++) + { + if(!rx[i]) continue; + e_max = fmax(e_max, beta[i]); + } - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); + log_verbose(" e_max=%.6lf\n", e_max); + log_verbose(" e_min=%.6lf\n", e_min); - rval = LP_relax(lp); - abort_if(rval, "LP_relax failed"); + if(DOUBLE_leq(obj, e_min)) + { + if(DOUBLE_geq(obj, e_max)) + *epsilon = INFINITY; + else + *epsilon = obj; - //rval = LP_write(lp, "scale-to-ahull.lp"); - //abort_if(rval, "LP_write failed"); - //UTIL_pause(); + goto CLEANUP; + } + else + { + for(int i = 0; i < nrays; i++) + if(rx[i] && DOUBLE_eq(beta[i], e_min)) + t[i] = 1; + } + } CLEANUP: - if(rmatind) free(rmatind); - if(rmatval) free(rmatval); + log_verbose(" e=%.6lf\n", *epsilon); + if(t) free(t); return rval; } @@ -1023,108 +1202,6 @@ 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) -{ - 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; - } - - rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, - rmatval); - abort_if(rval, "LP_add_rows failed"); - } - - // create constraint \sum_{r \in R} s_r = 1 - sense = 'E'; - rhs = 1.0; - - for(int i = 0; i < nrays; i++) - { - rmatind[i] = nrays + i; - rmatval[i] = 1.0; - } - - rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); - - // 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; -} - int GREEDY_ND_find_tight_rays(int nrows, int nrays, const double *f, @@ -1179,71 +1256,6 @@ 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) -{ - 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; -} - int GREEDY_ND_find_violated_cone(int nrows, int nrays, const double *f, @@ -1325,3 +1337,5 @@ CLEANUP: LP_free(&lp); return rval; } + +#endif // TEST_SOURCE \ No newline at end of file diff --git a/infinity/library/src/infinity-2d.c b/infinity/library/src/infinity-2d.c index b2e27c1..fff9315 100644 --- a/infinity/library/src/infinity-2d.c +++ b/infinity/library/src/infinity-2d.c @@ -524,29 +524,29 @@ CLEANUP: #ifndef TEST_SOURCE -int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) +int INFINITY_2D_generate_lfree(const struct MultiRowModel *model, + struct ConvLFreeSet *lfree) { - log_verbose("INFINITY_2D_generate_cut\n"); int rval = 0; int count = 0; int nrays = model->rays.nrays; - double *f = model->f; double *scale = 0; - double *rays = 0; - double lb[2], ub[2]; - for (int i = 0; i < nrays; i++) - bounds[i] = GREEDY_BIG_E; + double *f = lfree->f; + double *beta = lfree->beta; + double *rays = lfree->rays.values; + + lfree->nrows = 2; + lfree->rays.nrays = nrays; + memcpy(f, model->f, 2 * sizeof(double)); + memcpy(rays, model->rays.values, 2 * nrays * sizeof(double)); + for (int i = 0; i < nrays; i++) beta[i] = GREEDY_BIG_E; scale = (double*) malloc(nrays * sizeof(double)); - rays = (double*) malloc(2 * nrays * sizeof(double)); - abort_if(!rays, "could not allocate rays"); abort_if(!scale, "could not allocate scale"); - memcpy(rays, model->rays.values, 2 * nrays * sizeof(double)); - rval = scale_to_chull(rays, nrays, scale); abort_if(rval, "scale_to_chull failed"); @@ -558,7 +558,7 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) abort_if(count++ > 2 * nrays, "infinite loop"); - rval = get_bounding_box(2, nrays, rays, bounds, GREEDY_BIG_E, lb, ub); + rval = get_bounding_box(2, nrays, rays, beta, GREEDY_BIG_E, lb, ub); abort_if(rval, "get_bounding_box failed"); log_verbose(" box=[%.2lf %.2lf] [%.2lf %.2lf]\n", lb[0], ub[0], lb[1], ub[1]); @@ -581,7 +581,7 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) log_verbose(" p=%.2lf %.2lf\n", p[0], p[1]); - rval = bound(rays, bounds, nrays, f, p, &epsilon, v1, v2, &i1, &i2); + rval = bound(rays, beta, nrays, f, p, &epsilon, v1, v2, &i1, &i2); abort_if(rval, "bound failed"); log_verbose(" epsilon=%.2lf\n", epsilon); @@ -590,7 +590,7 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) { log_verbose(" found smaller epsilon: %.8lf\n", epsilon); - rval = get_bounding_box(2, nrays, rays, bounds, epsilon, lb, ub); + rval = get_bounding_box(2, nrays, rays, beta, epsilon, lb, ub); abort_if(rval, "get_bounding_box failed"); log_verbose(" p=%.2lf %.2lf\n", p[0], p[1]); @@ -625,11 +625,11 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) break; } - log_verbose(" updating bounds\n"); + log_verbose(" updating beta\n"); if(isinf(best_v1[0])) { - bounds[best_i1] = best_epsilon; - log_verbose(" bounds[%d]=%.8lf (exact)\n", best_i1, best_epsilon); + beta[best_i1] = best_epsilon; + log_verbose(" beta[%d]=%.8lf (exact)\n", best_i1, best_epsilon); } else { @@ -646,15 +646,15 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) if(!DOUBLE_geq(lambda, 0)) continue; - bounds[i] = fmin(bounds[i], lambda); - log_verbose(" bounds[%d]=%.8lf\n", i, bounds[i]); + beta[i] = fmin(beta[i], lambda); + log_verbose(" beta[%d]=%.8lf\n", i, beta[i]); } } //if(count > 0) //{ // for (int i = 0; i < nrays; i++) - // bounds[i] = fmin(bounds[i], best_epsilon); + // beta[i] = fmin(beta[i], best_epsilon); // break; //} @@ -663,7 +663,7 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) for (int k = 0; k < nrays; k++) { - if(bounds[k] < 100) continue; + if(beta[k] < 100) continue; is_split = 1; double *split_direction = &rays[2 * k]; @@ -700,13 +700,13 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) { const double *r = &rays[2 * i]; - lhs = (f[0] + r[0] * bounds[i]) * pi[0]; - lhs += (f[1] + r[1] * bounds[i]) * pi[1]; + lhs = (f[0] + r[0] * beta[i]) * pi[0]; + lhs += (f[1] + r[1] * beta[i]) * pi[1]; if (!(DOUBLE_leq(pi_zero, lhs) && DOUBLE_leq(lhs, pi_zero+1))) { log_verbose(" point %.4lf %.4lf falls outside of the split\n", - f[0] + r[0]*bounds[i], f[1] + r[1] * bounds[i]); + f[0] + r[0]*beta[i], f[1] + r[1] * beta[i]); is_split = 0; } } @@ -722,11 +722,10 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) } for(int i=0; invars; 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, model->rays.nrays, model->f, - model->rays.values, beta, &lp); + rval = GREEDY_create_psi_lp(nrows, lfree->rays.nrays, lfree->f, + lfree->rays.values, lfree->beta, &lp); abort_if(rval, "create_psi_lp failed"); + cut->nz = nvars; for(int i = 0; i < nvars; i++) { double value; @@ -117,15 +112,15 @@ static int create_cut(const struct Tableau *tableau, if(ENABLE_LIFTING && tableau->column_types[map->indices[i]] == MILP_INTEGER) { - rval = GREEDY_ND_pi(nrows, model->rays.nrays, model->f, - model->rays.values, beta, q, map->ray_scale[i], &lp, + rval = GREEDY_ND_pi(nrows, lfree->rays.nrays, lfree->f, + lfree->rays.values, lfree->beta, q, map->ray_scale[i], &lp, &value); abort_if(rval, "GREEDY_ND_pi failed"); } else { - rval = GREEDY_ND_psi(nrows, model->rays.nrays, model->f, - model->rays.values, beta, q, map->ray_scale[i], &lp, + rval = GREEDY_ND_psi(nrows, lfree->rays.nrays, lfree->f, + lfree->rays.values, lfree->beta, q, map->ray_scale[i], &lp, &value); abort_if(rval, "GREEDY_ND_psi failed"); } @@ -146,7 +141,7 @@ CLEANUP: return rval; } -static int select_rays(const struct RayMap *map, +static int filter_rays(const struct RayMap *map, const struct Tableau *tableau, struct MultiRowModel *model) { @@ -166,8 +161,8 @@ static int select_rays(const struct RayMap *map, for(int j = 0; j < (rays->nrays); j++) { double *q = LFREE_get_ray(&map->rays, j); - double norm = 0; + double norm = 0; for(int k = 0; k < nrows; k++) norm += fabs(r[k] - q[k]); @@ -182,7 +177,7 @@ static int select_rays(const struct RayMap *map, } log_debug(" norm_cutoff=%8.2lf nrays=%8d\n", norm_cutoff, - model->nrays); + model->rays.nrays); if(rays->nrays < MAX_N_RAYS) break; } @@ -253,7 +248,7 @@ CLEANUP: } static int write_sage_file(const struct MultiRowModel *model, - const double *beta, + const struct ConvLFreeSet *lfree, const char *filename) { int rval = 0; @@ -280,7 +275,7 @@ static int write_sage_file(const struct MultiRowModel *model, fprintf(fsage, "pi=vector([\n"); for(int k = 0; k < model->rays.nrays; k++) - fprintf(fsage, " %.12lf,\n", 1 / beta[k]); + fprintf(fsage, " %.12lf,\n", 1 / lfree->beta[k]); fprintf(fsage, "])\n"); CLEANUP: @@ -288,7 +283,8 @@ CLEANUP: return rval; } -static int dump_cut(const struct MultiRowModel *model, const double *beta) +static int dump_cut(const struct MultiRowModel *model, + const struct ConvLFreeSet *lfree) { int rval = 0; @@ -296,7 +292,7 @@ static int dump_cut(const struct MultiRowModel *model, const double *beta) sprintf(filename, "cut-%03d.sage", DUMP_CUT_N++); time_printf("Writing %s...\n", filename); - rval = write_sage_file(model, beta, filename); + rval = write_sage_file(model, lfree, filename); abort_if(rval, "write_sage_file failed"); CLEANUP: @@ -315,8 +311,8 @@ static int extract_model_from_tableau(const struct Tableau *tableau, 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"); + rval = filter_rays(map, tableau, model); + abort_if(rval, "filter_rays failed"); if(ENABLE_LIFTING) { @@ -339,16 +335,16 @@ CLEANUP: int INFINITY_generate_cut(const struct Tableau *tableau, struct Row *cut) { int rval = 0; - int nrows = tableau->nrows; int max_nrays = CG_total_nz(tableau); - double *beta = 0; + struct RayMap map; struct MultiRowModel model; - rval = CG_malloc_model(&model, nrows, max_nrays + 100); - abort_if(rval, "CG_malloc_model failed"); + struct ConvLFreeSet lfree; - struct RayMap map; - rval = CG_init_ray_map(&map, max_nrays, nrows); + rval = CG_init_model(&model, tableau->nrows, max_nrays + 100); + abort_if(rval, "CG_init_model failed"); + + rval = CG_init_ray_map(&map, max_nrays, tableau->nrows); abort_if(rval, "CG_init_ray_map failed"); rval = extract_model_from_tableau(tableau, &model, &map); @@ -357,16 +353,16 @@ int INFINITY_generate_cut(const struct Tableau *tableau, struct Row *cut) if(model.rays.nrays < 3) { rval = ERR_NO_CUT; - cut->pi = 0; - cut->indices = 0; goto CLEANUP; } - beta = (double *) malloc(model.rays.nrays * sizeof(double)); - abort_if(!beta, "could not allocate beta"); + rval = LFREE_init_conv(&lfree, tableau->nrows, model.rays.nrays); + abort_if(rval, "LFREE_init_conv failed"); - if(nrows == 2) rval = INFINITY_2D_generate_cut(&model, beta); - else rval = GREEDY_ND_generate_cut(&model, beta); + if(tableau->nrows == 2) + rval = INFINITY_2D_generate_lfree(&model, &lfree); + else + rval = INFINITY_ND_generate_lfree(&model, &lfree); if(rval) { @@ -376,15 +372,15 @@ int INFINITY_generate_cut(const struct Tableau *tableau, struct Row *cut) if(SHOULD_DUMP_CUTS) { - rval = dump_cut(&model, beta); + rval = dump_cut(&model, &lfree); abort_if(rval, "dump_cut failed"); } - rval = create_cut(tableau, &map, &model, beta, cut); - abort_if(rval, "create_cut failed"); + rval = create_cut_from_lfree(tableau, &map, &lfree, cut); + abort_if(rval, "create_cut_from_lfree failed"); CLEANUP: - if(beta) free(beta); + LFREE_free_conv(&lfree); CG_free_model(&model); CG_free_ray_map(&map); return rval; diff --git a/infinity/library/tests/greedy-nd-test.cpp b/infinity/library/tests/greedy-nd-test.cpp index 9ef10fa..1bbe643 100644 --- a/infinity/library/tests/greedy-nd-test.cpp +++ b/infinity/library/tests/greedy-nd-test.cpp @@ -25,6 +25,7 @@ extern "C" { #include #include #include +#include "../src/greedy-nd.c" } int ENABLE_LIFTING = 0; @@ -242,8 +243,8 @@ TEST(GreedyNDTest, bound_test_1) double epsilon; int tx[6]; - rval = GREEDY_ND_bound(2, 6, f, rays, x, beta1, &epsilon, tx); - abort_if(rval, "GREEDY_ND_bound failed"); + rval = bound(2, 6, f, rays, x, beta1, &epsilon, tx); + abort_if(rval, "bound failed"); EXPECT_NEAR(epsilon, 0.5, 1e-6); EXPECT_TRUE(tx[0]); EXPECT_FALSE(tx[1]); @@ -252,8 +253,8 @@ TEST(GreedyNDTest, bound_test_1) EXPECT_FALSE(tx[4]); EXPECT_FALSE(tx[5]); - rval = GREEDY_ND_bound(2, 6, f, rays, x, beta2, &epsilon, tx); - abort_if(rval, "GREEDY_ND_bound failed"); + rval = bound(2, 6, f, rays, x, beta2, &epsilon, tx); + abort_if(rval, "bound failed"); EXPECT_NEAR(epsilon, 1.0, 1e-6); EXPECT_TRUE(tx[0]); EXPECT_FALSE(tx[1]); @@ -262,8 +263,8 @@ TEST(GreedyNDTest, bound_test_1) EXPECT_TRUE(tx[4]); EXPECT_TRUE(tx[5]); - rval = GREEDY_ND_bound(2, 6, f, rays, x, beta3, &epsilon, tx); - abort_if(rval, "GREEDY_ND_bound failed"); + rval = bound(2, 6, f, rays, x, beta3, &epsilon, tx); + abort_if(rval, "bound failed"); EXPECT_EQ(epsilon, INFINITY); EXPECT_FALSE(tx[0]); EXPECT_FALSE(tx[1]); @@ -371,10 +372,9 @@ TEST(GreedyNDTest, generate_cut_test_1) double r3[] = { -1.0, 1.0 }; double r4[] = { 0.0, 1.0 }; double r5[] = { 1.0, 0.0 }; - double beta[6]; struct MultiRowModel model; - CG_malloc_model(&model, 2, 6); + CG_init_model(&model, 2, 6); LFREE_push_ray(&model.rays, r0); LFREE_push_ray(&model.rays, r1); LFREE_push_ray(&model.rays, r2); @@ -384,15 +384,18 @@ TEST(GreedyNDTest, generate_cut_test_1) model.f[0] = 0.5; model.f[1] = 0.5; - rval = GREEDY_ND_generate_cut(&model, beta); - abort_if(rval, "GREEDY_ND_generate_cut failed"); + struct ConvLFreeSet lfree; + LFREE_init_conv(&lfree, 2, 6); - EXPECT_NEAR(beta[0], 0.5, 1e-6); - EXPECT_NEAR(beta[1], 0.5, 1e-6); - EXPECT_NEAR(beta[2], 0.5, 1e-6); - EXPECT_NEAR(beta[3], 0.5, 1e-6); - EXPECT_NEAR(beta[4], 1.0, 1e-6); - EXPECT_NEAR(beta[5], 1.0, 1e-6); + rval = INFINITY_ND_generate_lfree(&model, &lfree); + abort_if(rval, "INFINITY_ND_generate_lfree failed"); + + EXPECT_NEAR(lfree.beta[0], 0.5, 1e-6); + EXPECT_NEAR(lfree.beta[1], 0.5, 1e-6); + EXPECT_NEAR(lfree.beta[2], 0.5, 1e-6); + EXPECT_NEAR(lfree.beta[3], 0.5, 1e-6); + EXPECT_NEAR(lfree.beta[4], 1.0, 1e-6); + EXPECT_NEAR(lfree.beta[5], 1.0, 1e-6); CLEANUP: if(rval) FAIL(); @@ -408,10 +411,9 @@ TEST(GreedyNDTest, generate_cut_test_2) double r3[] = { 0.0, -1.0, 0.0 }; double r4[] = { 0.0, 0.0, 1.0 }; double r5[] = { 0.0, 0.0, -1.0 }; - double beta[6]; struct MultiRowModel model; - CG_malloc_model(&model, 3, 6); + CG_init_model(&model, 3, 6); LFREE_push_ray(&model.rays, r0); LFREE_push_ray(&model.rays, r1); LFREE_push_ray(&model.rays, r2); @@ -422,15 +424,18 @@ TEST(GreedyNDTest, generate_cut_test_2) model.f[1] = 0.75; model.f[2] = 0.75; - rval = GREEDY_ND_generate_cut(&model, beta); - abort_if(rval, "GREEDY_ND_generate_cut failed"); + struct ConvLFreeSet lfree; + LFREE_init_conv(&lfree, 3, 6); + + rval = INFINITY_ND_generate_lfree(&model, &lfree); + abort_if(rval, "INFINITY_ND_generate_lfree failed"); - EXPECT_NEAR(beta[0], 0.75, 1e-6); - EXPECT_NEAR(beta[1], 2.25, 1e-6); - EXPECT_NEAR(beta[2], 0.75, 1e-6); - EXPECT_NEAR(beta[3], 2.25, 1e-6); - EXPECT_NEAR(beta[4], 0.75, 1e-6); - EXPECT_NEAR(beta[5], 2.25, 1e-6); + EXPECT_NEAR(lfree.beta[0], 0.75, 1e-6); + EXPECT_NEAR(lfree.beta[1], 2.25, 1e-6); + EXPECT_NEAR(lfree.beta[2], 0.75, 1e-6); + EXPECT_NEAR(lfree.beta[3], 2.25, 1e-6); + EXPECT_NEAR(lfree.beta[4], 0.75, 1e-6); + EXPECT_NEAR(lfree.beta[5], 2.25, 1e-6); CLEANUP: CG_free_model(&model); diff --git a/infinity/library/tests/infinity-2d-test.cpp b/infinity/library/tests/infinity-2d-test.cpp index 3609f59..6d7a100 100644 --- a/infinity/library/tests/infinity-2d-test.cpp +++ b/infinity/library/tests/infinity-2d-test.cpp @@ -30,7 +30,6 @@ extern "C" { TEST(Infinity2DTest, test_generate_cut_1) { int rval = 0; - double bounds[100]; double f[] = {1 / 4.0, 3 / 4.0}; double rays[] = { @@ -43,14 +42,17 @@ TEST(Infinity2DTest, test_generate_cut_1) const struct MultiRowModel model = {f , rays, 5, 2}; - rval = INFINITY_2D_generate_cut(&model, bounds); - abort_if(rval, "INFINITY_2D_generate_cut failed"); + struct ConvLFreeSet lfree; + LFREE_init_conv(&lfree, 2, 5); - EXPECT_NEAR(23 / 50.0, bounds[0], BOUNDS_EPSILON); - EXPECT_NEAR(23 / 42.0, bounds[1], BOUNDS_EPSILON); - EXPECT_NEAR(9 / 11.0, bounds[2], BOUNDS_EPSILON); - EXPECT_NEAR(9 / 11.0, bounds[3], BOUNDS_EPSILON); - EXPECT_NEAR(23 / 50.0, bounds[4], BOUNDS_EPSILON); + rval = INFINITY_2D_generate_lfree(&model, &lfree); + abort_if(rval, "INFINITY_2D_generate_lfree failed"); + + EXPECT_NEAR(23 / 50.0, lfree.beta[0], BOUNDS_EPSILON); + EXPECT_NEAR(23 / 42.0, lfree.beta[1], BOUNDS_EPSILON); + EXPECT_NEAR(9 / 11.0, lfree.beta[2], BOUNDS_EPSILON); + EXPECT_NEAR(9 / 11.0, lfree.beta[3], BOUNDS_EPSILON); + EXPECT_NEAR(23 / 50.0, lfree.beta[4], BOUNDS_EPSILON); CLEANUP: if (rval) FAIL(); @@ -59,7 +61,6 @@ TEST(Infinity2DTest, test_generate_cut_1) TEST(Infinity2DTest, test_generate_cut_2) { int rval = 0; - double bounds[100]; double f[] = {1 / 2.0, 1 / 2.0}; double rays[] = { -1.0, -1.0, @@ -70,15 +71,17 @@ TEST(Infinity2DTest, test_generate_cut_2) }; const struct MultiRowModel model = {f , rays, 5, 2}; + struct ConvLFreeSet lfree; + LFREE_init_conv(&lfree, 2, 5); - rval = INFINITY_2D_generate_cut(&model, bounds); - abort_if(rval, "INFINITY_2D_generate_cut failed"); + rval = INFINITY_2D_generate_lfree(&model, &lfree); + abort_if(rval, "INFINITY_2D_generate_lfree failed"); - EXPECT_NEAR(0.5, bounds[0], BOUNDS_EPSILON); - EXPECT_NEAR(0.5, bounds[1], BOUNDS_EPSILON); - EXPECT_NEAR(0.5, bounds[2], BOUNDS_EPSILON); - EXPECT_EQ(GREEDY_BIG_E, bounds[3]); - EXPECT_NEAR(0.5, bounds[4], BOUNDS_EPSILON); + EXPECT_NEAR(0.5, lfree.beta[0], BOUNDS_EPSILON); + EXPECT_NEAR(0.5, lfree.beta[1], BOUNDS_EPSILON); + EXPECT_NEAR(0.5, lfree.beta[2], BOUNDS_EPSILON); + EXPECT_EQ(GREEDY_BIG_E, lfree.beta[3]); + EXPECT_NEAR(0.5, lfree.beta[4], BOUNDS_EPSILON); CLEANUP: if (rval) FAIL(); @@ -87,18 +90,19 @@ TEST(Infinity2DTest, test_generate_cut_2) 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}; const struct MultiRowModel model = {f , rays, 3, 2}; + struct ConvLFreeSet lfree; + LFREE_init_conv(&lfree, 2, 5); - rval = INFINITY_2D_generate_cut(&model, bounds); - abort_if(rval, "INFINITY_2D_generate_cut failed"); + rval = INFINITY_2D_generate_lfree(&model, &lfree); + abort_if(rval, "INFINITY_2D_generate_lfree failed"); - EXPECT_NEAR(5.0, bounds[0], BOUNDS_EPSILON); - EXPECT_NEAR(17.0, bounds[2], BOUNDS_EPSILON); - EXPECT_EQ(GREEDY_BIG_E, bounds[1]); + EXPECT_NEAR(5.0, lfree.beta[0], BOUNDS_EPSILON); + EXPECT_NEAR(17.0, lfree.beta[2], BOUNDS_EPSILON); + EXPECT_EQ(GREEDY_BIG_E, lfree.beta[1]); CLEANUP: if (rval) FAIL(); @@ -277,8 +281,8 @@ TEST(Infinity2DTest, 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 = INFINITY_2D_generate_cut(rays, 6, f, bounds); -// abort_if(rval, "INFINITY_2D_generate_cut failed"); +// rval = INFINITY_2D_generate_lfree(rays, 6, f, bounds); +// abort_if(rval, "INFINITY_2D_generate_lfree failed"); // // EXPECT_NEAR(20.0, bounds[0], BOUNDS_EPSILON); // EXPECT_NEAR(20.0, bounds[1], BOUNDS_EPSILON); @@ -303,8 +307,8 @@ TEST(Infinity2DTest, find_containing_cone_test_3) // 0.04545454545454545581, 0.00000000000000000000, // 0.04545454545454545581, 0.02631578947368420907}; // -// rval = INFINITY_2D_generate_cut(rays, 6, f, bounds); -// abort_if(rval, "INFINITY_2D_generate_cut failed"); +// rval = INFINITY_2D_generate_lfree(rays, 6, f, bounds); +// abort_if(rval, "INFINITY_2D_generate_lfree failed"); // // EXPECT_NEAR(20.0, bounds[0], BOUNDS_EPSILON); // EXPECT_NEAR(20.0, bounds[1], BOUNDS_EPSILON); diff --git a/multirow/include/multirow/cg.h b/multirow/include/multirow/cg.h index 6dee19c..a364627 100644 --- a/multirow/include/multirow/cg.h +++ b/multirow/include/multirow/cg.h @@ -107,7 +107,7 @@ void CG_free_ray_map(struct RayMap *map); int CG_extract_f_from_tableau(const struct Tableau *tableau, double *f); -int CG_malloc_model(struct MultiRowModel *model, int nrows, int rays_capacity); +int CG_init_model(struct MultiRowModel *model, int nrows, int rays_capacity); void CG_free_model(struct MultiRowModel *model); diff --git a/multirow/include/multirow/lfree2d.h b/multirow/include/multirow/lfree2d.h index b186f34..608cd23 100644 --- a/multirow/include/multirow/lfree2d.h +++ b/multirow/include/multirow/lfree2d.h @@ -40,6 +40,15 @@ struct RayList int dim; }; + +struct ConvLFreeSet +{ + double *f; + struct RayList rays; + int nrows; + double *beta; +}; + int LFREE_2D_init(struct LFreeSet2D *set, int n_vertices, int n_lattice_points, @@ -69,4 +78,8 @@ void LFREE_free_ray_list(struct RayList *list); int LFREE_init_ray_list(struct RayList *list, int dim, int capacity); +int LFREE_init_conv(struct ConvLFreeSet *lfree, int dim, int max_nrays); + +void LFREE_free_conv(struct ConvLFreeSet *lfree); + #endif //LFREE_2D_H diff --git a/multirow/include/multirow/params.h b/multirow/include/multirow/params.h index 55be1bb..8ded0bc 100644 --- a/multirow/include/multirow/params.h +++ b/multirow/include/multirow/params.h @@ -31,7 +31,7 @@ * LOG_LEVEL_WARNING * LOG_LEVEL_ERROR */ -#define LOG_LEVEL LOG_LEVEL_INFO +#define LOG_LEVEL LOG_LEVEL_DEBUG /* * Maximum bounding-box size for naive algorithm @@ -48,7 +48,7 @@ */ #define N_RAYS 100 -#define ONLY_CUT -1 +#define ONLY_CUT 242 /* * Time limit for the computation (user time, in seconds). diff --git a/multirow/src/cg.c b/multirow/src/cg.c index c32f06c..5eba504 100644 --- a/multirow/src/cg.c +++ b/multirow/src/cg.c @@ -575,7 +575,10 @@ void CG_free(struct CG *cg) if (cg->tableau_rows) { for (int i = 0; i < cg->nrows; i++) + { LP_free_row(cg->tableau_rows[i]); + free(cg->tableau_rows[i]); + } free(cg->tableau_rows); } @@ -598,19 +601,16 @@ int CG_add_single_row_cuts(struct CG *cg, SingleRowGeneratorCallback generate) log_verbose("Generating cut %d...\n", i); - struct Row *cut = 0; - - cut = (struct Row *) malloc(sizeof(struct Row)); - abort_if(!cut, "could not allocate cut"); + struct Row cut; - rval = generate(row, cg->column_types, cut); + rval = generate(row, cg->column_types, &cut); abort_if(rval, "generate failed"); int ignored; - rval = add_cut(cg, cut, &ignored); + rval = add_cut(cg, &cut, &ignored); abort_if(rval, "add_cut failed"); - LP_free_row(cut); + LP_free_row(&cut); } CLEANUP: @@ -856,23 +856,24 @@ int CG_add_multirow_cuts(struct CG *cg, progress_increment(); } - struct Row *cut = 0; - - cut = (struct Row *) malloc(sizeof(struct Row)); - abort_if(!cut, "could not allocate cut"); + struct Tableau tableau = {nrows, rows, cg->column_types}; + struct Row cut; - cut->pi = 0; - cut->indices = 0; + int max_nz = CG_total_nz(&tableau); + cut.pi = (double *) malloc(max_nz * sizeof(double)); + cut.indices = (int *) malloc(max_nz * sizeof(int)); + abort_if(!cut.pi, "could not allocate cut.pi"); + abort_if(!cut.indices, "could not allocate cut.indices"); double initial_time = get_user_time(); - struct Tableau tableau = {nrows, rows, cg->column_types}; - rval = generate(&tableau, cut); + SHOULD_DUMP_CUTS = 1; + rval = generate(&tableau, &cut); if (rval == ERR_NO_CUT) { rval = 0; log_verbose("combination does not yield cut\n"); - LP_free_row(cut); + LP_free_row(&cut); goto NEXT_COMBINATION; } else abort_iff(rval, "generate failed (cut %d)", count); @@ -881,36 +882,36 @@ int CG_add_multirow_cuts(struct CG *cg, log_debug(" generate: %.2lf ms\n", elapsed_time * 1000); double dynamism; - rval = cut_dynamism(cut, &dynamism); + rval = cut_dynamism(&cut, &dynamism); abort_if(rval, "cut_dynamism failed"); if (dynamism > MAX_CUT_DYNAMISM) { log_debug("Discarding cut (dynamism=%.2lf)\n", dynamism); - LP_free_row(cut); + LP_free_row(&cut); goto NEXT_COMBINATION; } int ignored; - rval = add_cut(cg, cut, &ignored); + rval = add_cut(cg, &cut, &ignored); if (rval) { log_warn("invalid cut skipped (cut %d)\n", count); rval = 0; - LP_free_row(cut); + LP_free_row(&cut); goto NEXT_COMBINATION; } if_debug_level if (!ignored) - { - SHOULD_DUMP_CUTS = 1; - generate(&tableau, cut); - SHOULD_DUMP_CUTS = 0; - } + { + SHOULD_DUMP_CUTS = 1; + generate(&tableau, &cut); + SHOULD_DUMP_CUTS = 0; + } - LP_free_row(cut); + LP_free_row(&cut); } NEXT_COMBINATION: @@ -1019,7 +1020,7 @@ int CG_extract_f_from_tableau(const struct Tableau *tableau, double *f) return 0; } -int CG_malloc_model(struct MultiRowModel *model, int nrows, int rays_capacity) +int CG_init_model(struct MultiRowModel *model, int nrows, int rays_capacity) { int rval = 0; model->nrows = nrows; diff --git a/multirow/src/lfree2d.c b/multirow/src/lfree2d.c index b98e0b6..698a730 100644 --- a/multirow/src/lfree2d.c +++ b/multirow/src/lfree2d.c @@ -540,4 +540,29 @@ void LFREE_push_ray(struct RayList *list, const double *ray) double *dest = LFREE_get_ray(list, list->nrays); memcpy(dest, ray, list->dim * sizeof(double)); list->nrays++; +} + +int LFREE_init_conv(struct ConvLFreeSet *lfree, int dim, int max_nrays) +{ + int rval = 0; + + lfree->nrows = dim; + lfree->f = (double*) malloc(dim * sizeof(double)); + lfree->beta = (double*) malloc(max_nrays * sizeof(double)); + abort_if(!lfree->f, "could not allocate lfree->f"); + abort_if(!lfree->beta, "could not allocate lfree->beta"); + + rval = LFREE_init_ray_list(&lfree->rays, dim, max_nrays); + abort_if(rval, "LFREE_init_ray_list failed"); + +CLEANUP: + return rval; +} + +void LFREE_free_conv(struct ConvLFreeSet *lfree) +{ + if(!lfree) return; + free(lfree->f); + free(lfree->beta); + LFREE_free_ray_list(&lfree->rays); } \ No newline at end of file diff --git a/multirow/src/lp.c b/multirow/src/lp.c index 0372c49..b2f5496 100644 --- a/multirow/src/lp.c +++ b/multirow/src/lp.c @@ -80,10 +80,8 @@ void LP_free(struct LP *lp) void LP_free_row(struct Row *row) { if (!row) return; - if (row->pi) free(row->pi); if (row->indices) free(row->indices); - free(row); } int LP_add_rows(struct LP *lp, diff --git a/multirow/tests/cg-test.cpp b/multirow/tests/cg-test.cpp index 1e4a363..845367e 100644 --- a/multirow/tests/cg-test.cpp +++ b/multirow/tests/cg-test.cpp @@ -26,6 +26,7 @@ extern "C" { int BOOST_VAR = -1; double BOOST_FACTOR = 1.0; +int SHOULD_DUMP_CUTS = 0; TEST(CGTest, next_combination_test_1) { @@ -157,7 +158,7 @@ TEST(CGTest, extract_rays_from_rows_test) rval = CG_extract_rays_from_tableau(&tableau, &map); abort_if(rval, "CG_extract_rays_from_rows failed"); - EXPECT_EQ(rays.nrays, 4); + EXPECT_EQ(map.rays.nrays, 4); EXPECT_EQ(map.nvars, 7); EXPECT_DOUBLE_EQ(rays.values[0], 0.0);