From 208bc051bbba15e9c9a88460cc3ca9f29500b779 Mon Sep 17 00:00:00 2001 From: Alinson Xavier Date: Sun, 13 Aug 2017 13:25:13 -0400 Subject: [PATCH] infinity: Implement cone_bound_find_lambda --- infinity/library/src/infinity-nd.c | 36 +++++++ infinity/library/tests/infinity-nd-test.cpp | 112 ++++++++++++++++---- multirow/include/multirow/linalg.h | 8 +- multirow/include/multirow/util.h | 1 + multirow/src/linalg.c | 17 ++- multirow/tests/linalg-test.cpp | 26 ++++- 6 files changed, 163 insertions(+), 37 deletions(-) diff --git a/infinity/library/src/infinity-nd.c b/infinity/library/src/infinity-nd.c index ac2545b..8303d2c 100644 --- a/infinity/library/src/infinity-nd.c +++ b/infinity/library/src/infinity-nd.c @@ -22,6 +22,7 @@ #include #include +#include static long lp_count = 0; static double lp_time = 0; @@ -41,6 +42,41 @@ static double sfree_mip_time = 0; static long scale_ahull_lp_count = 0; static double scale_ahull_lp_time = 0; +static int cone_bound_find_lambda(const struct RayList *rays, + const double *f, + const double *x, + double *lambda) +{ + int rval = 0; + double *A = 0; + double *b = 0; + int dim = rays->dim; + int nrays = rays->nrays; + + A = (double*) malloc(nrays * dim * sizeof(double)); + b = (double*) malloc(dim * sizeof(double)); + abort_if(!A, "could not allocate A"); + abort_if(!b, "could not allocate b"); + + for (int i = 0; i < dim; i++) + { + b[i] = x[i] - f[i]; + for (int j = 0; j < nrays; j++) + A[i * nrays + j] = rays->values[j * dim + i]; + } + + rval = LINALG_solve(nrays, dim, A, b, lambda); + abort_if(rval, "LINALG_solve failed"); + + for (int i = 0; i < dim; i++) + abort_iff(lambda[0] < 0, "lambda[i] is negative (%lf)", i, lambda[i]); + +CLEANUP: + if(A) free(A); + if(b) free(b); + return rval; +} + static int create_sfree_mip(int nrows, int nrays, const double *f, diff --git a/infinity/library/tests/infinity-nd-test.cpp b/infinity/library/tests/infinity-nd-test.cpp index f3b1269..224d989 100644 --- a/infinity/library/tests/infinity-nd-test.cpp +++ b/infinity/library/tests/infinity-nd-test.cpp @@ -27,6 +27,7 @@ extern "C" { #include "../src/infinity-nd.c" } +const double E = 1e-6; TEST(InfinityNDTest, find_violated_cone_test) { @@ -164,11 +165,11 @@ TEST(InfinityNDTest, cone_bound_test_1) rval = cone_bound(2, 6, f, rays, rx1, x, beta, &epsilon); abort_if(rval, "cone_bound failed"); - EXPECT_NEAR(0.5, epsilon, 1e-6); + EXPECT_NEAR(0.5, epsilon, E); rval = cone_bound(2, 6, f, rays, rx2, x, beta, &epsilon); abort_if(rval, "cone_bound failed"); - EXPECT_NEAR(1.0, epsilon, 1e-6); + EXPECT_NEAR(1.0, epsilon, E); CLEANUP: if(rval) FAIL(); @@ -197,7 +198,7 @@ TEST(InfinityNDTest, cone_bound_test_2) rval = cone_bound(2, 2, f, rays, rx, x1, beta1, &epsilon); abort_if(rval, "cone_bound failed"); - EXPECT_NEAR(1.0, epsilon, 1e-6); + EXPECT_NEAR(1.0, epsilon, E); rval = cone_bound(2, 2, f, rays, rx, x1, beta2, &epsilon); abort_if(rval, "cone_bound failed"); @@ -205,7 +206,7 @@ TEST(InfinityNDTest, cone_bound_test_2) rval = cone_bound(2, 2, f, rays, rx, x2, beta2, &epsilon); abort_if(rval, "cone_bound failed"); - EXPECT_NEAR(1.0, epsilon, 1e-6); + EXPECT_NEAR(1.0, epsilon, E); rval = cone_bound(2, 2, f, rays, rx, x2, beta3, &epsilon); abort_if(rval, "cone_bound failed"); @@ -239,7 +240,7 @@ TEST(InfinityNDTest, bound_test_1) rval = bound(2, 6, f, rays, x, beta1, &epsilon, tx); abort_if(rval, "bound failed"); - EXPECT_NEAR(epsilon, 0.5, 1e-6); + EXPECT_NEAR(epsilon, 0.5, E); EXPECT_TRUE(tx[0]); EXPECT_FALSE(tx[1]); EXPECT_FALSE(tx[2]); @@ -249,7 +250,7 @@ TEST(InfinityNDTest, bound_test_1) rval = bound(2, 6, f, rays, x, beta2, &epsilon, tx); abort_if(rval, "bound failed"); - EXPECT_NEAR(epsilon, 1.0, 1e-6); + EXPECT_NEAR(epsilon, 1.0, E); EXPECT_TRUE(tx[0]); EXPECT_FALSE(tx[1]); EXPECT_FALSE(tx[2]); @@ -308,11 +309,11 @@ TEST(InfinityNDTest, psi_test) rval = INFINITY_psi(2, q1, 1.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); - EXPECT_NEAR(value, 2.0, 1e-6); + EXPECT_NEAR(value, 2.0, E); rval = INFINITY_psi(2, q2, 2.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); - EXPECT_NEAR(value, 8.0, 1e-6); + EXPECT_NEAR(value, 8.0, E); CLEANUP: LP_free(&lp); @@ -360,11 +361,11 @@ TEST(InfinityNDTest, psi_test_2) rval = INFINITY_psi(3, q1, 1.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); - EXPECT_NEAR(value, 1.0, 1e-6); + EXPECT_NEAR(value, 1.0, E); rval = INFINITY_psi(3, q2, 1.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); - EXPECT_NEAR(value, 2.0, 1e-6); + EXPECT_NEAR(value, 2.0, E); CLEANUP: if(rval) FAIL(); @@ -411,7 +412,7 @@ TEST(InfinityNDTest, psi_test_3) rval = INFINITY_psi(lfree.nrows, q, 1.0, &lp, &value); abort_if(rval, "GREDDY_ND_psi failed"); - EXPECT_NEAR(value, 1.0, 1e-6); + EXPECT_NEAR(value, 1.0, E); CLEANUP: if(rval) FAIL(); @@ -445,12 +446,12 @@ TEST(DISABLED_InfinityNDTest, generate_cut_test_1) 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); + EXPECT_NEAR(lfree.beta[0], 0.5, E); + EXPECT_NEAR(lfree.beta[1], 0.5, E); + EXPECT_NEAR(lfree.beta[2], 0.5, E); + EXPECT_NEAR(lfree.beta[3], 0.5, E); + EXPECT_NEAR(lfree.beta[4], 1.0, E); + EXPECT_NEAR(lfree.beta[5], 1.0, E); CLEANUP: if(rval) FAIL(); @@ -485,12 +486,12 @@ TEST(InfinityNDTest, generate_cut_test_2) rval = INFINITY_ND_generate_lfree(&model, &lfree); abort_if(rval, "INFINITY_ND_generate_lfree failed"); - 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); + EXPECT_NEAR(lfree.beta[0], 0.75, E); + EXPECT_NEAR(lfree.beta[1], 2.25, E); + EXPECT_NEAR(lfree.beta[2], 0.75, E); + EXPECT_NEAR(lfree.beta[3], 2.25, E); + EXPECT_NEAR(lfree.beta[4], 0.75, E); + EXPECT_NEAR(lfree.beta[5], 2.25, E); CLEANUP: CG_free_model(&model); @@ -534,3 +535,68 @@ TEST(InfinityNDTest, scale_to_ahull_test) CLEANUP: if(rval) FAIL(); } + +TEST(InfinityNDTest, cone_bound_find_lambda_test) +{ + int rval = 0; + + double f[] = { 1/2.0, 1/3.0 }; + double x[] = { 3.0, 3.0 }; + double ray_values[] = { + 1.0, 2.0, + 3.0, 2.0, + }; + + struct RayList rays = {.values = ray_values, .nrays = 2, .dim = 2}; + double lambda[2]; + + rval = cone_bound_find_lambda(&rays, f, x, lambda); + abort_if(rval, "cone_bound_find_lambda failed"); + + EXPECT_NEAR(lambda[0], 3/4.0, E); + EXPECT_NEAR(lambda[1], 7/12.0, E); + + CLEANUP: + if(rval) FAIL(); +} + +TEST(InfinityNDTest, cone_bound_find_lambda_test_2) +{ + int rval = 0; + + double f[] = { 1/2.0, 1/3.0, 5/6.0 }; + double x[] = { 3.0, 3.0, 6.0 }; + double ray_values[] = { + 1.0, 2.0, 3.0, + 3.0, 2.0, 5.0, + }; + + struct RayList rays = {.values = ray_values, .nrays = 2, .dim = 3}; + double lambda[2]; + + rval = cone_bound_find_lambda(&rays, f, x, lambda); + abort_if(rval, "cone_bound_find_lambda failed"); + + EXPECT_NEAR(lambda[0], 3/4.0, E); + EXPECT_NEAR(lambda[1], 7/12.0, E); + +CLEANUP: + if(rval) FAIL(); +} + +TEST(InfinityNDTest, cone_bound_find_lambda_test_3) +{ + int rval = 0; + double f[] = { 1/2.0, 1/2.0 }; + double x[] = { 3.0, 3.0 }; + double ray_values[] = { + 0.0, -1.0, + -1.0, 0.0, + }; + + struct RayList rays = {.values = ray_values, .nrays = 2, .dim = 2}; + double lambda[2]; + + rval = cone_bound_find_lambda(&rays, f, x, lambda); + EXPECT_NE(rval, 0); +} \ No newline at end of file diff --git a/multirow/include/multirow/linalg.h b/multirow/include/multirow/linalg.h index 142aca8..2e23961 100644 --- a/multirow/include/multirow/linalg.h +++ b/multirow/include/multirow/linalg.h @@ -33,8 +33,8 @@ double LINALG_dot(int n, double *x, double *y); double LINALG_norm(int n, double *x); /** - * Given an invertible n-by-n matrix A and an n-dimensional vector b, this - * function set x to A^(-1) b. Returns zero if the operation is successful and - * non-zero otherwise. + * Given a full rank m-by-n matrix A and an m-dimensional vector b, this + * function finds x such that Ax = b. Returns zero if the operation is + * successful and non-zero otherwise. */ -int LINALG_solve(int n, double *A, double *b, double *x); \ No newline at end of file +int LINALG_solve(int n, int m, double *A, double *b, double *x); \ No newline at end of file diff --git a/multirow/include/multirow/util.h b/multirow/include/multirow/util.h index e8156c9..fffca0c 100644 --- a/multirow/include/multirow/util.h +++ b/multirow/include/multirow/util.h @@ -78,6 +78,7 @@ void time_printf(const char *fmt, #define UNUSED(x) (void)(x) #define min(x,y) ((x) < (y) ? (x) : (y)) +#define max(x,y) ((x) > (y) ? (x) : (y)) #define sign(x) ((x) > 0 ? 1 : ((x) < 0 ? -1 : 0)) #define swap(x, y, T) do { T SWAP = x; x = y; y = SWAP; } while (0) diff --git a/multirow/src/linalg.c b/multirow/src/linalg.c index 21b34ee..67b3f48 100644 --- a/multirow/src/linalg.c +++ b/multirow/src/linalg.c @@ -17,8 +17,9 @@ * with this program. If not, see . */ -#include -#include +#include +#include +#include #include void LINALG_scale(int n, double *x, double alpha) @@ -36,16 +37,14 @@ double LINALG_norm(int n, double *x) return cblas_dasum(n, x, 1); } -int LINALG_solve(int n, double *A, double *b, double *x) +int LINALG_solve(int n, int m, double *A, double *b, double *x) { int rval = 0; - int ignored[n]; - double A_copy[n * n]; + double A_copy[n * m]; + memcpy(x, b, m * sizeof(double)); + memcpy(A_copy, A, n * m * sizeof(double)); - memcpy(x, b, n * sizeof(double)); - memcpy(A_copy, A, n * n * sizeof(double)); - - rval = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, A_copy, n, ignored, x, 1); + LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', m, n, 1, A_copy, n, x, 1); abort_if(rval, "LAPACKE_dgesv failed"); CLEANUP: diff --git a/multirow/tests/linalg-test.cpp b/multirow/tests/linalg-test.cpp index 2ddafec..a8e10a6 100644 --- a/multirow/tests/linalg-test.cpp +++ b/multirow/tests/linalg-test.cpp @@ -62,7 +62,7 @@ TEST(LinAlgTest, solve_test) double b[] = { 1.0, 3.0, 5.0 }; double x[] = { 0.0, 0.0, 0.0 }; - rval = LINALG_solve(3, A, b, x); + rval = LINALG_solve(3, 3, A, b, x); abort_if(rval, "LINALG_solve failed"); // Should compute x correctly @@ -78,6 +78,30 @@ TEST(LinAlgTest, solve_test) EXPECT_EQ(b[1], 3.0); EXPECT_EQ(b[2], 5.0); +CLEANUP: + if(rval) FAIL(); +} + +TEST(LinAlgTest, solve_test_2) +{ + int rval = 0; + + double A[] = { + 1.0, 3.0, + 0.0, 0.0, + 2.0, 2.0, + 3.0, 5.0, + }; + double b[] = { 5/2.0, 0.0, 8/3.0, 31/6.0 }; + double x[2]; + + rval = LINALG_solve(2, 4, A, b, x); + abort_if(rval, "LINALG_solve failed"); + + // Should compute x correctly + EXPECT_NEAR(x[0], 3/4.0, E); + EXPECT_NEAR(x[1], 7/12.0, E); + CLEANUP: if(rval) FAIL(); } \ No newline at end of file