From 3ded5cd96fa5ba2453c8de4035a4575824fe4c2b Mon Sep 17 00:00:00 2001 From: Alinson Xavier Date: Sat, 29 Apr 2017 23:49:45 -0400 Subject: [PATCH] Create structs (RayList, MultiRowModel, Tableau) --- infinity/library/include/infinity/greedy-nd.h | 22 +- infinity/library/include/infinity/infinity.h | 2 +- infinity/library/src/greedy-nd.c | 290 +++++++++-------- infinity/library/src/infinity-2d.c | 4 +- infinity/library/src/infinity.c | 292 ++++++------------ infinity/library/tests/greedy-nd-test.cpp | 62 ++-- infinity/library/tests/infinity-test.cpp | 45 ++- multirow/CMakeLists.txt | 1 + multirow/include/multirow/cg.h | 17 +- multirow/include/multirow/lfree2d.h | 22 +- multirow/src/cg.c | 53 ++-- multirow/src/lfree2d.c | 29 ++ multirow/tests/cg-test.cpp | 42 +-- 13 files changed, 425 insertions(+), 456 deletions(-) diff --git a/infinity/library/include/infinity/greedy-nd.h b/infinity/library/include/infinity/greedy-nd.h index ce295a3..f948f14 100644 --- a/infinity/library/include/infinity/greedy-nd.h +++ b/infinity/library/include/infinity/greedy-nd.h @@ -40,20 +40,16 @@ int GREEDY_ND_psi(const int nrows, double *value); 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); + 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 GREEDY_ND_generate_cut(int nrows, - int nrays, - const double *f, - const double *rays, - double *beta); +int GREEDY_ND_generate_cut(const struct MultiRowModel *model, double *beta); int GREEDY_ND_bound(int nrows, int nrays, diff --git a/infinity/library/include/infinity/infinity.h b/infinity/library/include/infinity/infinity.h index 1b08782..74e0803 100644 --- a/infinity/library/include/infinity/infinity.h +++ b/infinity/library/include/infinity/infinity.h @@ -19,6 +19,6 @@ #include #include -int INFINITY_generate_cut(struct Tableau *tableau, struct Row *cut); +int INFINITY_generate_cut(const struct Tableau *tableau, struct Row *cut); #endif //MULTIROW_INFINITY_H diff --git a/infinity/library/src/greedy-nd.c b/infinity/library/src/greedy-nd.c index 059ff42..bbca1d9 100644 --- a/infinity/library/src/greedy-nd.c +++ b/infinity/library/src/greedy-nd.c @@ -21,7 +21,6 @@ #include #include -#include #include static long lp_count = 0; @@ -56,7 +55,7 @@ int find_interior_point_enum(const int nrows, abort_if(nrows != 3, "not implemented"); double *beta2 = 0; - beta2 = (double*) malloc(nrays * sizeof(double)); + beta2 = (double *) malloc(nrays * sizeof(double)); abort_if(!beta2, "could not allocate beta2"); for(int i = 0; i < nrays; i++) @@ -65,7 +64,7 @@ int find_interior_point_enum(const int nrows, } struct LP lp; - + rval = LP_open(&lp); abort_if(rval, "LP_open failed"); @@ -81,9 +80,12 @@ int find_interior_point_enum(const int nrows, for(int x3 = -M; x3 <= M; x3++) { double value; - double q[3] = {x1 - f[0], x2 - f[1], x3 - f[2]}; + double q[3] = {x1 - f[0], + x2 - f[1], + x3 - f[2]}; - rval = GREEDY_ND_psi(nrows, nrays, f, rays, beta2, q, 1, &lp, &value); + rval = GREEDY_ND_psi(nrows, nrays, f, rays, beta2, q, 1, &lp, + &value); abort_if(rval, "GREEDY_ND_psi failed"); if(value < best_value) @@ -102,6 +104,7 @@ CLEANUP: LP_free(&lp); return rval; } + int find_interior_point_cplex(const int nrows, const int nrays, const double *f, @@ -113,7 +116,7 @@ int find_interior_point_cplex(const int nrows, { int rval = 0; struct LP lp; - + rval = LP_open(&lp); abort_if(rval, "LP_open failed"); @@ -145,11 +148,11 @@ int find_interior_point_cplex(const int nrows, if_verbose_level { for(int i = 0; i < nrows; i++) - log_verbose(" x%d = %.8lf\n", i, x[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]); + if(x[i + nrows] > 0.00001) + log_verbose(" t%d = %.8lf\n", i, x[i + nrows]); } rval = LP_get_obj_val(&lp, &objval); @@ -184,7 +187,7 @@ int GREEDY_create_psi_lp(const int nrows, int rval = 0; int rmatbeg = 0; - int* rmatind = 0; + int *rmatind = 0; double *rmatval = 0; rmatind = (int *) malloc(nrays * sizeof(int)); @@ -196,20 +199,20 @@ int GREEDY_create_psi_lp(const int nrows, abort_if(rval, "LP_create failed"); // create lambda variables - for (int i = 0; i < nrays; i++) + 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++) + for(int j = 0; j < nrows; j++) { char sense = 'E'; double rhs = 0; int nz = 0; - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) { const double *ri = &rays[i * nrows]; rmatind[nz] = i; @@ -217,8 +220,7 @@ int GREEDY_create_psi_lp(const int nrows, nz++; } - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); abort_if(rval, "LP_add_rows failed"); } @@ -232,20 +234,20 @@ int GREEDY_create_psi_lp(const int nrows, } CLEANUP: - if (rmatind) free(rmatind); - if (rmatval) free(rmatval); + 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) + 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; @@ -259,7 +261,8 @@ int GREEDY_ND_pi(const int nrows, 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 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, @@ -277,11 +280,9 @@ int GREEDY_ND_pi(const int nrows, 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 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, @@ -339,26 +340,24 @@ CLEANUP: return rval; } -int GREEDY_ND_generate_cut(int nrows, - int nrays, - const double *f, - const double *rays, - double *beta) +int GREEDY_ND_generate_cut(const struct MultiRowModel *model, double *beta) { - log_debug("GREEDY_ND_generate_cut\n"); int rval = 0; + int nrows = model->nrows; + int nrays = model->rays.nrays; + double *f = model->f; double *x = 0; int *t = 0; int *tx = 0; - - t = (int*) malloc(nrays * sizeof(int)); - tx = (int*) malloc(nrays * sizeof(int)); + + 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"); - x = (double*) malloc((nrows + nrays) * sizeof(double)); + x = (double *) malloc((nrows + nrays) * sizeof(double)); abort_if(!x, "could not allocate x"); for(int i = 0; i < nrays; i++) @@ -386,11 +385,10 @@ int GREEDY_ND_generate_cut(int nrows, long x_count = 0; double epsilon; - double previous_epsilon = -INFINITY; for(int i = 0; i < nrows; i++) - log_verbose(" f[%d] = %.12lf\n", i, f[i]); - + log_verbose(" f[%d] = %.12lf\n", i, f[i]); + while(1) { it++; @@ -413,13 +411,15 @@ int GREEDY_ND_generate_cut(int nrows, if(nrows == 3) { - rval = find_interior_point_enum(nrows, nrays, f, rays, beta, epsilon, x, &found); + rval = find_interior_point_enum(nrows, nrays, f, + model->rays.values, beta, epsilon, x, &found); abort_if(rval, "find_interior_point_enum failed"); } if(!found) { - rval = find_interior_point_cplex(nrows, nrays, f, rays, beta, epsilon, x, &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; @@ -432,8 +432,8 @@ int GREEDY_ND_generate_cut(int nrows, abort_if(x_count > 1000, "infinite loop"); double epsilon_x; - rval = GREEDY_ND_bound(nrows, nrays, f, rays, x, beta, &epsilon_x, - tx); + 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; @@ -469,10 +469,10 @@ int GREEDY_ND_generate_cut(int nrows, else if(!skip_ahull) { double alpha; - const double *d = &rays[i * nrows]; + const double *d = LFREE_get_ray(&model->rays, i); - rval = GREEDY_ND_scale_to_ahull(nrows, nrays, rays, t, beta, - epsilon, d, &alpha); + 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)) @@ -488,33 +488,38 @@ int GREEDY_ND_generate_cut(int nrows, log_debug(" beta[%2d] = %.4lf\n", i, beta[i]); } - previous_epsilon = epsilon; 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", + 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, + 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: @@ -540,9 +545,9 @@ int GREEDY_ND_bound(int nrows, double *fbar = 0; double *sbar = 0; - rx = (int*) malloc(nrays * sizeof(int)); - fbar = (double*) malloc(nrows * sizeof(double)); - sbar = (double*) malloc(nrays * sizeof(double)); + 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"); @@ -558,7 +563,7 @@ int GREEDY_ND_bound(int nrows, abort_if(count > 100, "infinite loop"); rval = GREEDY_ND_find_violated_cone(nrows, nrays, f, rays, x, beta, - *epsilon, rx, sbar, &found); + *epsilon, rx, sbar, &found); abort_if(rval, "GREEDY_ND_find_violated_cone failed"); if(!found) break; @@ -588,7 +593,6 @@ int GREEDY_ND_bound(int nrows, abort_if(prev_epsilon < *epsilon, "epsilon should never increase"); } - for(int i = 0; i < nrays; i++) tx[i] = 0; @@ -599,7 +603,8 @@ int GREEDY_ND_bound(int nrows, } else { - rval = GREEDY_ND_find_tight_rays(nrows, nrays, fbar, rays, x, beta, *epsilon, tx); + rval = GREEDY_ND_find_tight_rays(nrows, nrays, fbar, rays, x, beta, + *epsilon, tx); abort_if(rval, "GREEDY_ND_find_tight_rays failed"); } @@ -610,7 +615,6 @@ CLEANUP: return rval; } - static int create_find_epsilon_lp(int nrows, int nrays, const double *f, @@ -629,10 +633,10 @@ static int create_find_epsilon_lp(int nrows, int nz = 0; int *map = 0; int rmatbeg = 0; - int* rmatind = 0; + int *rmatind = 0; double *rmatval = 0; - map = (int*) malloc(nrays * sizeof(int)); + map = (int *) malloc(nrays * sizeof(int)); abort_if(!map, "could not allocate map"); rmatind = (int *) malloc((nrays + 1 + nrows) * sizeof(int)); @@ -646,7 +650,7 @@ static int create_find_epsilon_lp(int nrows, // create lambda variables int rx_count = 0; - for (int i = 0; i < nrays + 1; i++) + for(int i = 0; i < nrays + 1; i++) { if(i < nrays && !rx[i]) continue; @@ -669,20 +673,20 @@ static int create_find_epsilon_lp(int nrows, log_verbose("rx_count=%d\n", rx_count); // create y variables - for (int i = 0; i < nrows; i++) + 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"); } // create constraint y = \lambda_x x + \sum_{t \in T} \lambda_r (f + \beta_r r) - for (int j = 0; j < nrows; j++) + for(int j = 0; j < nrows; j++) { sense = 'E'; rhs = 0.0; nz = 0; - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) { if(!t[i]) continue; const double *ri = &rays[i * nrows]; @@ -700,25 +704,24 @@ static int create_find_epsilon_lp(int nrows, rmatval[nz] = -1.0; nz++; - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + 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++) + for(int j = 0; j < nrows; j++) { sense = 'E'; rhs = f[j]; nz = 0; - for (int i = 0; i < nrays; i++) + 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]; + rmatval[nz] = -ri[j]; nz++; } @@ -726,8 +729,7 @@ static int create_find_epsilon_lp(int nrows, rmatval[nz] = 1.0; nz++; - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); abort_if(rval, "LP_add_rows failed"); } @@ -736,7 +738,7 @@ static int create_find_epsilon_lp(int nrows, rhs = 1.0; nz = 0; - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) { if(!t[i]) continue; rmatind[nz] = map[i]; @@ -748,8 +750,7 @@ static int create_find_epsilon_lp(int nrows, rmatval[nz] = 1.0; nz++; - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); abort_if(rval, "LP_add_rows failed"); rval = LP_relax(lp); @@ -761,8 +762,8 @@ static int create_find_epsilon_lp(int nrows, CLEANUP: if(map) free(map); - if (rmatind) free(rmatind); - if (rmatval) free(rmatval); + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); return rval; } @@ -781,7 +782,7 @@ int GREEDY_ND_cone_bound(int nrows, int *t = 0; long it = 0; - t = (int*) malloc(nrays * sizeof(int)); + t = (int *) malloc(nrays * sizeof(int)); abort_if(!t, "could not allocate t"); for(int i = 0; i < nrays; i++) @@ -798,7 +799,8 @@ int GREEDY_ND_cone_bound(int nrows, rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = create_find_epsilon_lp(nrows, nrays, f, rays, t, rx, x, beta, &lp); + rval = create_find_epsilon_lp(nrows, nrays, f, rays, t, rx, x, beta, + &lp); abort_if(rval, "create_find_epsilon_lp failed"); int infeasible; @@ -874,7 +876,6 @@ CLEANUP: return rval; } - static int create_scale_to_ahull_lp(int nrows, int nrays, const double *rays, @@ -891,7 +892,7 @@ static int create_scale_to_ahull_lp(int nrows, int nz; int rmatbeg = 0; - int* rmatind = 0; + int *rmatind = 0; double *rmatval = 0; rmatind = (int *) malloc((nrays + 1) * sizeof(int)); @@ -908,14 +909,14 @@ static int create_scale_to_ahull_lp(int nrows, // create lambda variables - for (int i = 0; i < nrays; i++) + 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++) + for(int j = 0; j < nrows; j++) { sense = 'E'; rhs = 0.0; @@ -925,20 +926,19 @@ static int create_scale_to_ahull_lp(int nrows, rmatval[nz] = d[j]; nz++; - for (int i = 0; i < nrays; i++) + 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]; + 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); + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); abort_if(rval, "LP_add_rows failed"); } @@ -947,7 +947,7 @@ static int create_scale_to_ahull_lp(int nrows, rhs = 1.0; nz = 0; - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) { if(!rx[i]) continue; @@ -956,8 +956,7 @@ static int create_scale_to_ahull_lp(int nrows, nz++; } - rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + rval = LP_add_rows(lp, 1, nz, &rhs, &sense, &rmatbeg, rmatind, rmatval); abort_if(rval, "LP_add_rows failed"); rval = LP_relax(lp); @@ -968,8 +967,8 @@ static int create_scale_to_ahull_lp(int nrows, //UTIL_pause(); CLEANUP: - if (rmatind) free(rmatind); - if (rmatval) free(rmatval); + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); return rval; } @@ -989,15 +988,15 @@ int GREEDY_ND_scale_to_ahull(int nrows, struct LP lp; double *x = 0; - x = (double*) malloc((nrays + 1) * sizeof(double)); + x = (double *) malloc((nrays + 1) * sizeof(double)); abort_if(!x, "could not allocate x"); double initial_time = get_user_time(); rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = create_scale_to_ahull_lp(nrows, nrays, rays, rx, beta, epsilon, - d, &lp); + rval = create_scale_to_ahull_lp(nrows, nrays, rays, rx, beta, epsilon, d, + &lp); abort_if(rval, "create_scale_to_ahull_lp failed"); int infeasible; @@ -1024,7 +1023,6 @@ CLEANUP: return rval; } - static int create_tight_rays_lp(int nrows, int nrays, const double *f, @@ -1041,7 +1039,7 @@ static int create_tight_rays_lp(int nrows, char sense; int rmatbeg = 0; - int* rmatind = 0; + int *rmatind = 0; double *rmatval = 0; rmatind = (int *) malloc(nrays * sizeof(int)); @@ -1053,26 +1051,26 @@ static int create_tight_rays_lp(int nrows, abort_if(rval, "LP_create failed"); // create lambda variables - for (int i = 0; i < nrays; i++) + 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++) + 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++) + for(int j = 0; j < nrows; j++) { sense = 'E'; rhs = x[j] - f[j]; - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) { const double *ri = &rays[i * nrows]; rmatind[i] = nrays + i; @@ -1081,7 +1079,7 @@ static int create_tight_rays_lp(int nrows, } rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + rmatval); abort_if(rval, "LP_add_rows failed"); } @@ -1089,18 +1087,17 @@ static int create_tight_rays_lp(int nrows, sense = 'E'; rhs = 1.0; - for (int i = 0; i < nrays; i++) + 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); + 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++) + for(int i = 0; i < nrays; i++) { sense = 'G'; rhs = delta; @@ -1111,8 +1108,7 @@ static int create_tight_rays_lp(int nrows, rmatind[1] = nrays + i; rmatval[1] = 1.0; - rval = LP_add_rows(lp, 1, 2, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + rval = LP_add_rows(lp, 1, 2, &rhs, &sense, &rmatbeg, rmatind, rmatval); abort_if(rval, "LP_add_rows failed"); } @@ -1124,8 +1120,8 @@ static int create_tight_rays_lp(int nrows, CLEANUP: - if (rmatind) free(rmatind); - if (rmatval) free(rmatval); + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); return rval; } @@ -1145,15 +1141,15 @@ int GREEDY_ND_find_tight_rays(int nrows, struct LP lp; double *sbar = 0; - sbar = (double*) malloc(2 * nrays * sizeof(double)); + sbar = (double *) malloc(2 * nrays * sizeof(double)); abort_if(!sbar, "could not allocate sbar"); double initial_time = get_user_time(); rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = create_tight_rays_lp(nrows, nrays, f, rays, x, beta, epsilon, - delta, &lp); + rval = create_tight_rays_lp(nrows, nrays, f, rays, x, beta, epsilon, delta, + &lp); abort_if(rval, "create_tight_rays_lp failed"); int infeasible = 0; @@ -1171,11 +1167,11 @@ int GREEDY_ND_find_tight_rays(int nrows, rval = LP_get_x(&lp, sbar); abort_if(rval, "LP_get_x failed"); - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) tx[i] = DOUBLE_iszero(sbar[i]); - for (int i = 0; i < nrays; i++) - log_verbose(" tx[%d]=%d\n", i, tx[i]); + for(int i = 0; i < nrays; i++) + log_verbose(" tx[%d]=%d\n", i, tx[i]); CLEANUP: if(sbar) free(sbar); @@ -1183,7 +1179,6 @@ CLEANUP: return rval; } - static int create_violated_cone_lp(int nrows, int nrays, const double *f, @@ -1199,7 +1194,7 @@ static int create_violated_cone_lp(int nrows, char sense; int rmatbeg = 0; - int* rmatind = 0; + int *rmatind = 0; double *rmatval = 0; rmatind = (int *) malloc(nrays * sizeof(int)); @@ -1211,19 +1206,19 @@ static int create_violated_cone_lp(int nrows, abort_if(rval, "LP_create failed"); // create s variables - for (int i = 0; i < nrays; i++) + 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++) + for(int j = 0; j < nrows; j++) { sense = 'E'; rhs = x[j] - f[j]; - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) { const double *ri = &rays[i * nrows]; rmatind[i] = i; @@ -1232,7 +1227,7 @@ static int create_violated_cone_lp(int nrows, } rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, - rmatval); + rmatval); abort_if(rval, "LP_add_rows failed"); } @@ -1244,8 +1239,8 @@ static int create_violated_cone_lp(int nrows, //UTIL_pause(); CLEANUP: - if (rmatind) free(rmatind); - if (rmatval) free(rmatval); + if(rmatind) free(rmatind); + if(rmatval) free(rmatval); return rval; } @@ -1270,7 +1265,8 @@ int GREEDY_ND_find_violated_cone(int nrows, rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = create_violated_cone_lp(nrows, nrays, f, rays, x, beta, epsilon, &lp); + rval = create_violated_cone_lp(nrows, nrays, f, rays, x, beta, epsilon, + &lp); abort_if(rval, "create_violated_cone_lp failed"); int infeasible; @@ -1286,7 +1282,7 @@ int GREEDY_ND_find_violated_cone(int nrows, rval = LP_get_x(&lp, sbar); abort_if(rval, "LP_get_x failed"); - for(int i=0; i 1e-9); @@ -1320,7 +1316,7 @@ int GREEDY_ND_find_violated_cone(int nrows, 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]); + log_verbose(" r[%d]=%.8lf %.8lf\n", i, m * r[0], m * r[1]); } } } diff --git a/infinity/library/src/infinity-2d.c b/infinity/library/src/infinity-2d.c index af3d2f8..b2e27c1 100644 --- a/infinity/library/src/infinity-2d.c +++ b/infinity/library/src/infinity-2d.c @@ -529,7 +529,7 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) log_verbose("INFINITY_2D_generate_cut\n"); int rval = 0; int count = 0; - int nrays = model->nrays; + int nrays = model->rays.nrays; double *f = model->f; double *scale = 0; @@ -545,7 +545,7 @@ int INFINITY_2D_generate_cut(const struct MultiRowModel *model, double *bounds) abort_if(!rays, "could not allocate rays"); abort_if(!scale, "could not allocate scale"); - memcpy(rays, model->rays, 2 * nrays * sizeof(double)); + memcpy(rays, model->rays.values, 2 * nrays * sizeof(double)); rval = scale_to_chull(rays, nrays, scale); abort_if(rval, "scale_to_chull failed"); diff --git a/infinity/library/src/infinity.c b/infinity/library/src/infinity.c index 68e0bd6..038d7eb 100644 --- a/infinity/library/src/infinity.c +++ b/infinity/library/src/infinity.c @@ -22,8 +22,8 @@ #include #include -#include #include +#include struct SortPair { @@ -51,112 +51,43 @@ static int _qsort_cmp_rays_angle(const void *p1, const void *p2) /** * Sorts a list of rays according to their angle. - * - * @param rays the rays to be sorted - * @param nrays the number of rays in the list - * @param beta a list of doubles, to be sorted together with the rays - * - * @return zero if successful, non-zero otherwise */ -static int sort_rays_angle(double *rays, int nrays, double *beta) +static int sort_rays_by_angle(struct RayList *rays) { int rval = 0; - + int nrays = rays->nrays; double *rays_copy = 0; - double *beta_copy = 0; struct SortPair *pairs = 0; pairs = (struct SortPair *) malloc(nrays * sizeof(struct SortPair)); rays_copy = (double *) malloc(2 * nrays * sizeof(double)); - beta_copy = (double *) malloc(nrays * sizeof(double)); abort_if(!pairs, "could not allocate pairs"); abort_if(!rays_copy, "could not allocate rays_copy"); - abort_if(!beta_copy, "could not allocate beta_copy"); - memcpy(rays_copy, rays, 2 * nrays * sizeof(double)); - memcpy(beta_copy, beta, nrays * sizeof(double)); + memcpy(rays_copy, rays->values, 2 * nrays * sizeof(double)); - for (int i = 0; i < nrays; i++) + for(int i = 0; i < nrays; i++) { pairs[i].index = i; - pairs[i].data = &rays[2 * i]; + pairs[i].data = &rays->values[2 * i]; } qsort(pairs, (size_t) nrays, sizeof(struct SortPair), _qsort_cmp_rays_angle); - for (int i = 0; i < nrays; i++) - { - beta[i] = beta_copy[pairs[i].index]; - memcpy(&rays[2 * i], &rays_copy[2 * pairs[i].index], - 2 * sizeof(double)); - } + rays->nrays = 0; + for(int i = 0; i < nrays; i++) + LFREE_push_ray(rays, &rays_copy[2 * pairs[i].index]); CLEANUP: - if (pairs) free(pairs); - if (rays_copy) free(rays_copy); - if (beta_copy) free(beta_copy); + free(pairs); + free(rays_copy); return rval; } -static void print_row(const struct Row *row) -{ - time_printf("Row:\n"); - for (int i = 0; i < row->nz; i++) - time_printf(" %.4lfx%d\n", row->pi[i], row->indices[i]); - time_printf(" <= %.4lf [%d]\n", row->pi_zero, row->head); -} - -static void print_rays(const struct Tableau *tableau, - const struct RayMap *map, - const double *f, - const double *rays, - int nrays) -{ - int nrows = tableau->nrows; - - time_printf("Ray map:\n"); - for (int i = 0; i < map->nvars; i++) - time_printf(" %4d: %4d x%-4d (scale=%.4lf)\n", i, - map->variable_to_ray[i], map->indices[i], map->ray_scale[i]); - - time_printf("Origin:\n"); - for (int i = 0; i < nrows; i++) - time_printf(" %20.12lf\n", f[i]); - - time_printf("Rays:\n"); - for (int i = 0; i < nrays; i++) - { - time_printf(" "); - - for (int j = 0; j < nrows; j++) - printf("%20.12lf ", rays[i * nrows + j]); - - printf("angle=%.4lf ", atan2(rays[i * nrows], rays[i * nrows + 1])); - printf("norm=%.4lf ", - fabs(rays[i * nrows]) + fabs(rays[i * nrows + 1])); - - printf("[ "); - for (int j = 0; j < map->nvars; j++) - if (map->variable_to_ray[j] == i) - printf("%d ", map->indices[j]); - printf("]\n"); - } -} - -static void print_cut(const struct Row *cut) -{ - time_printf("Generated cut:\n"); - for (int i = 0; i < cut->nz; i++) - time_printf(" %.4lfx%d\n", cut->pi[i], cut->indices[i]); - time_printf(" <= %.4lf\n", cut->pi_zero); -} - static int create_cut(const struct Tableau *tableau, const struct RayMap *map, - const double *f, - const int lfree_nrays, - const double *lfree_rays, + const struct MultiRowModel *model, const double *beta, struct Row *cut) { @@ -171,29 +102,31 @@ static int create_cut(const struct Tableau *tableau, abort_if(!cut->indices, "could not allocate cut->indices"); struct LP lp; - rval = LP_open(&lp); abort_if(rval, "LP_open failed"); - rval = GREEDY_create_psi_lp(nrows, lfree_nrays, f, lfree_rays, beta, &lp); + rval = GREEDY_create_psi_lp(nrows, model->rays.nrays, model->f, + model->rays.values, beta, &lp); abort_if(rval, "create_psi_lp failed"); - for (int i = 0; i < nvars; i++) + for(int i = 0; i < nvars; i++) { double value; - const double *q = &map->rays[map->variable_to_ray[i] * nrows]; + const double *q = LFREE_get_ray(&map->rays, map->variable_to_ray[i]); - if (ENABLE_LIFTING && + if(ENABLE_LIFTING && tableau->column_types[map->indices[i]] == MILP_INTEGER) { - rval = GREEDY_ND_pi(nrows, lfree_nrays, f, lfree_rays, beta, q, - map->ray_scale[i], &lp, &value); + rval = GREEDY_ND_pi(nrows, model->rays.nrays, model->f, + model->rays.values, beta, q, map->ray_scale[i], &lp, + &value); abort_if(rval, "GREEDY_ND_pi failed"); } else { - rval = GREEDY_ND_psi(nrows, lfree_nrays, f, lfree_rays, beta, q, - map->ray_scale[i], &lp, &value); + rval = GREEDY_ND_psi(nrows, model->rays.nrays, model->f, + model->rays.values, beta, q, map->ray_scale[i], &lp, + &value); abort_if(rval, "GREEDY_ND_psi failed"); } @@ -213,57 +146,53 @@ CLEANUP: return rval; } -static int select_rays(const struct RayMap map, +static int select_rays(const struct RayMap *map, const struct Tableau *tableau, struct MultiRowModel *model) { int rval = 0; int nrows = tableau->nrows; + struct RayList *rays = &model->rays; - for (double norm_cutoff = 0.00; norm_cutoff <= 5.0; norm_cutoff += 0.1) + for(double norm_cutoff = 0.00; norm_cutoff <= 5.0; norm_cutoff += 0.1) { - model->nrays = 0; + rays->nrays = 0; - for (int i = 0; i < map.nrays; i++) + for(int i = 0; i < map->rays.nrays; i++) { int keep = 1; - double *r = &map.rays[tableau->nrows * i]; + double *r = LFREE_get_ray(&map->rays, i); - for (int j = 0; j < (model->nrays); j++) + for(int j = 0; j < (rays->nrays); j++) { - double *q = &map.rays[tableau->nrows * j]; + double *q = LFREE_get_ray(&map->rays, j); double norm = 0; - for (int k = 0; k < nrows; k++) + for(int k = 0; k < nrows; k++) norm += fabs(r[k] - q[k]); - if (norm <= norm_cutoff) + if(norm <= norm_cutoff) { keep = 0; break; } } - if (keep) - { - memcpy(&model->rays[nrows * (model->nrays)], r, - nrows * sizeof(double)); - model->nrays++; - } + if(keep) LFREE_push_ray(rays, r); } log_debug(" norm_cutoff=%8.2lf nrays=%8d\n", norm_cutoff, model->nrays); - if (model->nrays < MAX_N_RAYS) break; + if(rays->nrays < MAX_N_RAYS) break; } CLEANUP: return rval; } -static int -append_extra_rays(const struct Tableau *tableau, double *rays, int *nrays) +static int append_extra_rays(const struct Tableau *tableau, + struct RayList *rays) { int rval = 0; int nrows = tableau->nrows; @@ -272,7 +201,7 @@ append_extra_rays(const struct Tableau *tableau, double *rays, int *nrays) abort_if(nrows > 3, "not implemented"); - if (nrows == 2) + if(nrows == 2) { extra_rays[0] = 0.0001; extra_rays[1] = 0.0000; @@ -286,7 +215,7 @@ append_extra_rays(const struct Tableau *tableau, double *rays, int *nrays) n_extra_rays = 3; } - if (nrows == 3) + if(nrows == 3) { extra_rays[0] = 0.0000; extra_rays[1] = 0.0000; @@ -307,58 +236,50 @@ append_extra_rays(const struct Tableau *tableau, double *rays, int *nrays) n_extra_rays = 4; } - for (int i = 0; i < n_extra_rays; i++) + for(int i = 0; i < n_extra_rays; i++) { double *r = &extra_rays[nrows * i]; double scale; int found, index; - rval = CG_find_ray(nrows, rays, *nrays, r, &found, &scale, &index); + rval = CG_find_ray(rays, r, &found, &scale, &index); abort_if(rval, "CG_find_ray failed"); - if (!found) - { - memcpy(&rays[nrows * (*nrays)], r, nrows * - sizeof(double)); - (*nrays)++; - } + if(!found) LFREE_push_ray(rays, r); } CLEANUP: return rval; } -static int write_sage_file(int nrows, - int nrays, - const double *f, - const double *rays, +static int write_sage_file(const struct MultiRowModel *model, const double *beta, const char *filename) { int rval = 0; + int nrows = model->nrows; FILE *fsage = fopen(filename, "w"); abort_iff(!fsage, "could not open %s", filename); fprintf(fsage, "f=vector(["); - for (int i = 0; i < nrows; i++) - fprintf(fsage, "%.20lf,", f[i]); + for(int i = 0; i < nrows; i++) + fprintf(fsage, "%.20lf,", model->f[i]); fprintf(fsage, "])\n"); fprintf(fsage, "R=matrix([\n"); - for (int i = 0; i < nrays; i++) + for(int i = 0; i < model->rays.nrays; i++) { + double *r = LFREE_get_ray(&model->rays, i); fprintf(fsage, " ["); - for (int j = 0; j < nrows; j++) - { - fprintf(fsage, "%.20lf,", rays[i * nrows + j]); - } + for(int j = 0; j < nrows; j++) + fprintf(fsage, "%.20lf,", r[j]); fprintf(fsage, "],\n"); } fprintf(fsage, "])\n"); fprintf(fsage, "pi=vector([\n"); - for (int k = 0; k < nrays; k++) + for(int k = 0; k < model->rays.nrays; k++) fprintf(fsage, " %.12lf,\n", 1 / beta[k]); fprintf(fsage, "])\n"); @@ -367,11 +288,7 @@ CLEANUP: return rval; } -static int dump_cut(const struct Tableau *tableau, - const double *rays, - int nrays, - const double *f, - const double *beta) +static int dump_cut(const struct MultiRowModel *model, const double *beta) { int rval = 0; @@ -379,52 +296,65 @@ static int dump_cut(const struct Tableau *tableau, sprintf(filename, "cut-%03d.sage", DUMP_CUT_N++); time_printf("Writing %s...\n", filename); - rval = write_sage_file(tableau->nrows, nrays, f, rays, beta, filename); + rval = write_sage_file(model, beta, filename); abort_if(rval, "write_sage_file failed"); CLEANUP: return rval; } -#ifndef TEST_SOURCE - -int INFINITY_generate_cut(struct Tableau *tableau, struct Row *cut) +static int extract_model_from_tableau(const struct Tableau *tableau, + struct MultiRowModel *model, + struct RayMap *map) { int rval = 0; - int max_nrays = 0; - int nrows = tableau->nrows; - double *beta = 0; - - for (int i = 0; i < tableau->nrows; i++) - max_nrays += tableau->rows[i]->nz; - struct MultiRowModel model; - rval = CG_init_model(&model, nrows, max_nrays + 100); - abort_if(rval, "CG_init_model failed"); - - rval = CG_extract_f_from_tableau(tableau, model.f); + rval = CG_extract_f_from_tableau(tableau, model->f); abort_if(rval, "CG_extract_f_from_tableau failed"); - struct RayMap map; - rval = CG_init_ray_map(&map, max_nrays, nrows); - abort_if(rval, "CG_init_ray_map failed"); - - rval = CG_extract_rays_from_tableau(tableau, &map); + rval = CG_extract_rays_from_tableau(tableau, map); abort_if(rval, "CG_extract_rays_from_rows failed"); - rval = select_rays(map, tableau, &model); + rval = select_rays(map, tableau, model); abort_if(rval, "select_rays failed"); - if (ENABLE_LIFTING) + if(ENABLE_LIFTING) { - rval = append_extra_rays(tableau, model.rays, &model.nrays); + rval = append_extra_rays(tableau, &model->rays); abort_if(rval, "append_extra_rays failed"); } - beta = (double *) malloc(model.nrays * sizeof(double)); - abort_if(!beta, "could not allocate beta"); + if(tableau->nrows == 2) + { + rval = sort_rays_by_angle(&model->rays); + abort_if(rval, "sort_rays_by_angle failed"); + } - if (model.nrays < 3) +CLEANUP: + return rval; +} + +#ifndef TEST_SOURCE + +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 MultiRowModel model; + rval = CG_malloc_model(&model, nrows, max_nrays + 100); + abort_if(rval, "CG_malloc_model failed"); + + struct RayMap map; + rval = CG_init_ray_map(&map, max_nrays, nrows); + abort_if(rval, "CG_init_ray_map failed"); + + rval = extract_model_from_tableau(tableau, &model, &map); + abort_if(rval, "extract_model_from_tableau failed"); + + if(model.rays.nrays < 3) { rval = ERR_NO_CUT; cut->pi = 0; @@ -432,47 +362,29 @@ int INFINITY_generate_cut(struct Tableau *tableau, struct Row *cut) goto CLEANUP; } - log_debug("Selected %d rays\n", nrays); - if_verbose_level print_rays(tableau, &map, model.f, map.rays, map.nrays); + beta = (double *) malloc(model.rays.nrays * sizeof(double)); + abort_if(!beta, "could not allocate beta"); - log_verbose("Computing lattice-free set...\n"); - if (nrows == 2) - { - rval = sort_rays_angle(model.rays, model.nrays, beta); - abort_if(rval, "sort_rays_angle failed"); + if(nrows == 2) rval = INFINITY_2D_generate_cut(&model, beta); + else rval = GREEDY_ND_generate_cut(&model, beta); - rval = INFINITY_2D_generate_cut(&model, beta); - if (rval) - { - rval = ERR_NO_CUT; - goto CLEANUP; - } - } - else + if(rval) { - rval = GREEDY_ND_generate_cut(nrows, model.nrays, model.f, model.rays, - beta); - if (rval) - { - rval = ERR_NO_CUT; - goto CLEANUP; - } + rval = ERR_NO_CUT; + goto CLEANUP; } - if (SHOULD_DUMP_CUTS) + if(SHOULD_DUMP_CUTS) { - rval = dump_cut(tableau, model.rays, model.nrays, model.f, beta); + rval = dump_cut(&model, beta); abort_if(rval, "dump_cut failed"); } - rval = create_cut(tableau, &map, model.f, model.nrays, model.rays, beta, - cut); + rval = create_cut(tableau, &map, &model, beta, cut); abort_if(rval, "create_cut failed"); - if_verbose_level print_cut(cut); - CLEANUP: - if (beta) free(beta); + if(beta) free(beta); 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 fea9360..9ef10fa 100644 --- a/infinity/library/tests/greedy-nd-test.cpp +++ b/infinity/library/tests/greedy-nd-test.cpp @@ -22,6 +22,8 @@ extern "C" { #include #include #include +#include +#include #include } @@ -363,19 +365,26 @@ TEST(GreedyNDTest, generate_cut_test_1) { int rval = 0; - double f[] = { 0.5, 0.5 }; - double rays[] = - { - 1.0, 1.0, - 1.0, -1.0, - -1.0, -1.0, - -1.0, 1.0, - 0.0, 1.0, - 1.0, 0.0 - }; + double r0[] = { 1.0, 1.0 }; + double r1[] = { 1.0, -1.0 }; + double r2[] = { -1.0, -1.0 }; + double r3[] = { -1.0, 1.0 }; + double r4[] = { 0.0, 1.0 }; + double r5[] = { 1.0, 0.0 }; double beta[6]; - rval = GREEDY_ND_generate_cut(2, 6, f, rays, beta); + struct MultiRowModel model; + CG_malloc_model(&model, 2, 6); + LFREE_push_ray(&model.rays, r0); + LFREE_push_ray(&model.rays, r1); + LFREE_push_ray(&model.rays, r2); + LFREE_push_ray(&model.rays, r3); + LFREE_push_ray(&model.rays, r4); + LFREE_push_ray(&model.rays, r5); + 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"); EXPECT_NEAR(beta[0], 0.5, 1e-6); @@ -393,19 +402,27 @@ TEST(GreedyNDTest, generate_cut_test_2) { int rval = 0; - double f[] = { 0.75, 0.75, 0.75}; - double rays[] = - { - 1.0, 0.0, 0.0, - -1.0, 0.0, 0.0, - 0.0, 1.0, 0.0, - 0.0, -1.0, 0.0, - 0.0, 0.0, 1.0, - 0.0, 0.0, -1.0 - }; + double r0[] = { 1.0, 0.0, 0.0 }; + double r1[] = {-1.0, 0.0, 0.0 }; + double r2[] = { 0.0, 1.0, 0.0 }; + 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]; - rval = GREEDY_ND_generate_cut(3, 6, f, rays, beta); + struct MultiRowModel model; + CG_malloc_model(&model, 3, 6); + LFREE_push_ray(&model.rays, r0); + LFREE_push_ray(&model.rays, r1); + LFREE_push_ray(&model.rays, r2); + LFREE_push_ray(&model.rays, r3); + LFREE_push_ray(&model.rays, r4); + LFREE_push_ray(&model.rays, r5); + model.f[0] = 0.75; + 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"); EXPECT_NEAR(beta[0], 0.75, 1e-6); @@ -416,6 +433,7 @@ TEST(GreedyNDTest, generate_cut_test_2) EXPECT_NEAR(beta[5], 2.25, 1e-6); CLEANUP: + CG_free_model(&model); if(rval) FAIL(); } diff --git a/infinity/library/tests/infinity-test.cpp b/infinity/library/tests/infinity-test.cpp index 812009a..3a0dd79 100644 --- a/infinity/library/tests/infinity-test.cpp +++ b/infinity/library/tests/infinity-test.cpp @@ -53,39 +53,36 @@ TEST(InfinityTest, sort_rays_angle_test) { int rval = 0; - int n_rays = 5; - double beta[] = {0, 1, 2, 3, 4}; - SortPair sp0, sp1, sp2, sp3, sp4; + double r0[] = { 1.0, 1.0 }; + double r1[] = { 1.0, 0.0 }; + double r2[] = { 1.0, -1.0 }; + double r3[] = { -1.0, 0.0 }; + double r4[] = { 2.0, 0.0 }; - double rays[] = { - 1.0, 1.0, - 1.0, 0.0, - 1.0, -1.0, - -1.0, 0.0, - 2.0, 0.0, - }; + RayList rays; + LFREE_init_ray_list(&rays, 2, 5); + LFREE_push_ray(&rays, r0); + LFREE_push_ray(&rays, r1); + LFREE_push_ray(&rays, r2); + LFREE_push_ray(&rays, r3); + LFREE_push_ray(&rays, r4); - rval = sort_rays_angle(rays, n_rays, beta); - abort_if(rval, "sort_rays_angle failed"); + rval = sort_rays_by_angle(&rays); + abort_if(rval, "sort_rays_by_angle failed"); - sp0 = { 0, &rays[0] }; - sp1 = { 1, &rays[2] }; - sp2 = { 2, &rays[4] }; - sp3 = { 3, &rays[6] }; - sp4 = { 4, &rays[8] }; + SortPair sp0, sp1, sp2, sp3, sp4; + sp0 = { 0, LFREE_get_ray(&rays, 0) }; + sp1 = { 1, LFREE_get_ray(&rays, 1) }; + sp2 = { 2, LFREE_get_ray(&rays, 2) }; + sp3 = { 3, LFREE_get_ray(&rays, 3) }; + sp4 = { 4, LFREE_get_ray(&rays, 4) }; EXPECT_LE(_qsort_cmp_rays_angle(&sp0, &sp1), 0); EXPECT_LE(_qsort_cmp_rays_angle(&sp1, &sp2), 0); EXPECT_LE(_qsort_cmp_rays_angle(&sp2, &sp3), 0); EXPECT_LE(_qsort_cmp_rays_angle(&sp3, &sp4), 0); - EXPECT_EQ(beta[0], 3); - EXPECT_EQ(beta[1], 0); - EXPECT_TRUE(beta[2] == 1 || beta[2] == 4); - EXPECT_TRUE(beta[3] == 1 || beta[3] == 4); - EXPECT_TRUE(beta[2] != beta[3]); - EXPECT_EQ(beta[4], 2); - CLEANUP: + LFREE_free_ray_list(&rays); if (rval) FAIL(); } diff --git a/multirow/CMakeLists.txt b/multirow/CMakeLists.txt index ab7d940..ff7a311 100644 --- a/multirow/CMakeLists.txt +++ b/multirow/CMakeLists.txt @@ -13,6 +13,7 @@ set(COMMON_SOURCES include/multirow/geometry.h include/multirow/lfree2d.h include/multirow/lp.h + include/multirow/mir.h include/multirow/rational.h include/multirow/stats.h include/multirow/params.h diff --git a/multirow/include/multirow/cg.h b/multirow/include/multirow/cg.h index e373077..6dee19c 100644 --- a/multirow/include/multirow/cg.h +++ b/multirow/include/multirow/cg.h @@ -17,7 +17,8 @@ #ifndef MULTIROW_CG_H #define MULTIROW_CG_H -#include "lp.h" +#include +#include struct CG { @@ -47,8 +48,7 @@ struct Tableau struct RayMap { - int nrays; - double *rays; + struct RayList rays; int *variable_to_ray; double *ray_scale; int *indices; @@ -58,8 +58,7 @@ struct RayMap struct MultiRowModel { double *f; - double *rays; - int nrays; + struct RayList rays; int nrows; }; @@ -96,9 +95,7 @@ int CG_boost_variable(int var, int *indices, int nz); -int CG_find_ray(int dim, - const double *rays, - int nrays, +int CG_find_ray(const struct RayList *rays, const double *r, int *found, double *scale, @@ -110,8 +107,10 @@ void CG_free_ray_map(struct RayMap *map); int CG_extract_f_from_tableau(const struct Tableau *tableau, double *f); -int CG_init_model(struct MultiRowModel *model, int nrows, int max_nrays); +int CG_malloc_model(struct MultiRowModel *model, int nrows, int rays_capacity); void CG_free_model(struct MultiRowModel *model); +int CG_total_nz(const struct Tableau *tableau); + #endif //MULTIROW_CG_H diff --git a/multirow/include/multirow/lfree2d.h b/multirow/include/multirow/lfree2d.h index f1f6278..b186f34 100644 --- a/multirow/include/multirow/lfree2d.h +++ b/multirow/include/multirow/lfree2d.h @@ -1,4 +1,3 @@ - /* Copyright (c) 2015 Alinson Xavier * * This program is free software: you can redistribute it and/or modify @@ -18,6 +17,8 @@ #ifndef LFREE_2D_H #define LFREE_2D_H +#include + struct LFreeSet2D { double f[2]; @@ -32,6 +33,13 @@ struct LFreeSet2D double *halfspaces; }; +struct RayList +{ + double *values; + int nrays; + int dim; +}; + int LFREE_2D_init(struct LFreeSet2D *set, int n_vertices, int n_lattice_points, @@ -51,8 +59,14 @@ int LFREE_2D_print_set(const struct LFreeSet2D *set); int LFREE_2D_translate_set(struct LFreeSet2D *set, double dx, double dy); -int LFREE_2D_get_bounding_box(const struct LFreeSet2D *set, - int *lb, - int *ub); +int LFREE_2D_get_bounding_box(const struct LFreeSet2D *set, int *lb, int *ub); + +void LFREE_push_ray(struct RayList *list, const double *ray); + +double* LFREE_get_ray(const struct RayList *list, int index); + +void LFREE_free_ray_list(struct RayList *list); + +int LFREE_init_ray_list(struct RayList *list, int dim, int capacity); #endif //LFREE_2D_H diff --git a/multirow/src/cg.c b/multirow/src/cg.c index 6911c56..c32f06c 100644 --- a/multirow/src/cg.c +++ b/multirow/src/cg.c @@ -371,9 +371,7 @@ CLEANUP: * parallel to the given one. Returns whether a matching ray was found, * the index of the matching ray and the scaling factor. */ -int CG_find_ray(int dim, - const double *rays, - int nrays, +int CG_find_ray(const struct RayList *rays, const double *r, int *found, double *scale, @@ -381,10 +379,12 @@ int CG_find_ray(int dim, { *found = 0; - for (int i = 0; i < nrays; i++) + for (int i = 0; i < rays->nrays; i++) { + double *q = LFREE_get_ray(rays, i); + int match; - check_rays_parallel(dim, r, &rays[dim * i], &match, scale); + check_rays_parallel(rays->dim, r, q, &match, scale); if (match) { @@ -404,11 +404,9 @@ int CG_extract_rays_from_tableau(const struct Tableau *tableau, int rval = 0; int nrows = tableau->nrows; - double *rays = map->rays; struct Row **rows = tableau->rows; map->nvars = 0; - map->nrays = 0; int *i = 0; int *idx = 0; @@ -423,7 +421,7 @@ int CG_extract_rays_from_tableau(const struct Tableau *tableau, while (1) { - double *r = &rays[nrows * (map->nrays)]; + double *r = LFREE_get_ray(&map->rays, map->rays.nrays); int idx_min = INT_MAX; @@ -462,20 +460,19 @@ int CG_extract_rays_from_tableau(const struct Tableau *tableau, for (int j = 0; j < nrows; j++) log_verbose(" r[%d] = %.12lf\n", j, r[j]); - rval = CG_find_ray(nrows, rays, map->nrays, r, &found, &scale, - &ray_index); + rval = CG_find_ray(&map->rays, r, &found, &scale, &ray_index); abort_if(rval, "CG_find_ray failed"); if (!found) { log_verbose(" ray is new\n"); scale = 1.0; - ray_index = (map->nrays)++; + ray_index = map->rays.nrays++; } else { log_verbose(" ray equals:\n"); - double *q = &rays[ray_index * nrows]; + double *q = LFREE_get_ray(&map->rays, ray_index); for (int j = 0; j < nrows; j++) log_verbose(" r[%d] = %.12lf\n", j, q[j]); } @@ -488,10 +485,9 @@ int CG_extract_rays_from_tableau(const struct Tableau *tableau, NEXT_RAY:; } - for (int j = 0; j < map->nrays; j++) + for (int j = 0; j < map->rays.nrays; j++) { - double *r = &rays[nrows * j]; - + double *r = LFREE_get_ray(&map->rays, j); double max_scale = 0.0; for (int k = 0; k < map->nvars; k++) @@ -992,11 +988,12 @@ int CG_init_ray_map(struct RayMap *map, int max_nrays, int nrows) map->variable_to_ray = (int *) malloc(max_nrays * sizeof(int)); map->indices = (int *) malloc(max_nrays * sizeof(int)); map->ray_scale = (double *) malloc(max_nrays * sizeof(double)); - map->rays = (double *) malloc(nrows * max_nrays * sizeof(double)); abort_if(!map->variable_to_ray, "could not allocate variable_to_ray"); abort_if(!map->indices, "could not allocate indices"); abort_if(!map->ray_scale, "could not allocate ray_scale"); - abort_if(!map->rays, "could not allocate rays"); + + rval = LFREE_init_ray_list(&map->rays, nrows, max_nrays); + abort_if(rval, "LFREE_init_ray_list failed"); CLEANUP: return rval; @@ -1008,7 +1005,7 @@ void CG_free_ray_map(struct RayMap *map) free(map->variable_to_ray); free(map->indices); free(map->ray_scale); - free(map->rays); + LFREE_free_ray_list(&map->rays); } int CG_extract_f_from_tableau(const struct Tableau *tableau, double *f) @@ -1022,16 +1019,16 @@ int CG_extract_f_from_tableau(const struct Tableau *tableau, double *f) return 0; } -int CG_init_model(struct MultiRowModel *model, int nrows, int max_nrays) +int CG_malloc_model(struct MultiRowModel *model, int nrows, int rays_capacity) { int rval = 0; - - model->nrays = 0; model->nrows = nrows; + + rval = LFREE_init_ray_list(&model->rays, nrows, rays_capacity); + abort_if(rval, "LFREE_init_ray_list failed"); + model->f = (double*) malloc(nrows * sizeof(double)); - model->rays = (double*) malloc(max_nrays * sizeof(double)); abort_if(!model->f, "could not allocate f"); - abort_if(!model->rays, "could not allocate rays"); CLEANUP: return rval; @@ -1040,8 +1037,16 @@ CLEANUP: void CG_free_model(struct MultiRowModel *model) { if(!model) return; - free(model->rays); free(model->f); + LFREE_free_ray_list(&model->rays); +} + +int CG_total_nz(const struct Tableau *tableau) +{ + int total_nz = 0; + for(int i = 0; i < tableau->nrows; i++) + total_nz += tableau->rows[i]->nz; + return total_nz; } #endif // TEST_SOURCE diff --git a/multirow/src/lfree2d.c b/multirow/src/lfree2d.c index 490c6a0..b98e0b6 100644 --- a/multirow/src/lfree2d.c +++ b/multirow/src/lfree2d.c @@ -512,3 +512,32 @@ CLEANUP: return rval; } +int LFREE_init_ray_list(struct RayList *list, int dim, int capacity) +{ + int rval = 0; + list->nrays = 0; + list->dim = dim; + list->values = (double*) malloc(capacity * dim * sizeof(double)); + abort_if(!list->values, "could not allocate list->values"); + +CLEANUP: + return rval; +} + +void LFREE_free_ray_list(struct RayList *list) +{ + if(!list) return; + free(list->values); +} + +double* LFREE_get_ray(const struct RayList *list, int index) +{ + return &list->values[index * list->dim]; +} + +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++; +} \ No newline at end of file diff --git a/multirow/tests/cg-test.cpp b/multirow/tests/cg-test.cpp index 18d6524..1e4a363 100644 --- a/multirow/tests/cg-test.cpp +++ b/multirow/tests/cg-test.cpp @@ -71,7 +71,7 @@ TEST(CGTest, next_combination_test_2) EXPECT_EQ(count, 10); - CLEANUP: +CLEANUP: if(rval) FAIL(); } @@ -109,6 +109,7 @@ TEST(CGTest, extract_rays_from_rows_test) { int rval = 0; + char column_types[16] = {0}; double pi1[] = { 1.0, 1.0, 1.0, 2.0, 1.0 }; int indices1[] = { 1, 7, 8, 12, 14 }; struct Row row1 = @@ -143,36 +144,37 @@ TEST(CGTest, extract_rays_from_rows_test) }; struct Row *rows[] = { &row1, &row2, &row3 }; + struct Tableau tableau = {3, rows, column_types}; - int nz; - int nrays; int indices[1000]; int variable_to_ray[1000]; - double rays[1000]; double ray_scale[1000]; - rval = CG_extract_rays_from_tableau(3, rows, rays, &nrays, variable_to_ray, - ray_scale, indices, &nz); + struct RayList rays; + LFREE_init_ray_list(&rays, 3, 1000); + struct RayMap map = {rays, variable_to_ray, ray_scale, indices, 0}; + + rval = CG_extract_rays_from_tableau(&tableau, &map); abort_if(rval, "CG_extract_rays_from_rows failed"); - EXPECT_EQ(nrays, 4); - EXPECT_EQ(nz, 7); + EXPECT_EQ(rays.nrays, 4); + EXPECT_EQ(map.nvars, 7); - EXPECT_DOUBLE_EQ(rays[0], 0.0); - EXPECT_DOUBLE_EQ(rays[1], 0.0); - EXPECT_DOUBLE_EQ(rays[2], -1.0); + EXPECT_DOUBLE_EQ(rays.values[0], 0.0); + EXPECT_DOUBLE_EQ(rays.values[1], 0.0); + EXPECT_DOUBLE_EQ(rays.values[2], -1.0); - EXPECT_DOUBLE_EQ(rays[3], 0.0); - EXPECT_DOUBLE_EQ(rays[4], -1.0); - EXPECT_DOUBLE_EQ(rays[5], 0.0); + EXPECT_DOUBLE_EQ(rays.values[3], 0.0); + EXPECT_DOUBLE_EQ(rays.values[4], -1.0); + EXPECT_DOUBLE_EQ(rays.values[5], 0.0); - EXPECT_DOUBLE_EQ(rays[6], -2.0); - EXPECT_DOUBLE_EQ(rays[7], -2.0); - EXPECT_DOUBLE_EQ(rays[8], 0.0); + EXPECT_DOUBLE_EQ(rays.values[6], -2.0); + EXPECT_DOUBLE_EQ(rays.values[7], -2.0); + EXPECT_DOUBLE_EQ(rays.values[8], 0.0); - EXPECT_DOUBLE_EQ(rays[ 9], 0.0); - EXPECT_DOUBLE_EQ(rays[10], -1.0); - EXPECT_DOUBLE_EQ(rays[11], -1.0); + EXPECT_DOUBLE_EQ(rays.values[ 9], 0.0); + EXPECT_DOUBLE_EQ(rays.values[10], -1.0); + EXPECT_DOUBLE_EQ(rays.values[11], -1.0); EXPECT_EQ(indices[0], 3); EXPECT_EQ(indices[1], 5);