isoron
/
multirow
Archived
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

New project structure

Browse Source
selection
Alinson S. Xavier 9 years ago
commit 7098f65110
50 changed files with 11246 additions and 0 deletions
Show Stats Download Patch File Download Diff File

3
.gitmodules vendored
Unescape View File

@ -0,0 +1,3 @@
[submodule "googletest"]
path = googletest
url = https://github.com/google/googletest.git

20
CMakeLists.txt
Unescape View File

@ -0,0 +1,20 @@
cmake_minimum_required(VERSION 2.8)
project(multirow)
set(CMAKE_LIBRARY_PATH ${CMAKE_LIBRARY_PATH} /Users/axavier/Applications/IBM/ILOG/CPLEX_Studio1262/cplex/lib/x86-64_osx/static_pic)
list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/cmake")
find_package(CPLEX REQUIRED)
include_directories(${CPLEX_INCLUDE_DIR})
include_directories(multirow/include)
include_directories(lifting/library/include)
include_directories(infinity/library/include)
add_subdirectory(googletest/googletest)
include_directories(${gtest_SOURCE_DIR}/include)
add_subdirectory(multirow)
add_subdirectory(lifting/library)
add_subdirectory(lifting/benchmark)
add_subdirectory(infinity/library)
add_subdirectory(infinity/benchmark)

17
cmake/FindCPLEX.cmake
Unescape View File

@ -0,0 +1,17 @@
find_library(CPLEX_LIBRARIES
NAMES cplex cplex1220 cplex1240 cplex1260 cplex1261 cplex1262)
find_path(CPLEX_INCLUDE_DIR
NAMES ilcplex/cplex.h
PATHS /Users/axavier/Applications/IBM/ILOG/CPLEX_Studio1262/cplex/include)
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(
CPLEX
DEFAULT_MSG
CPLEX_LIBRARIES
CPLEX_INCLUDE_DIR)
mark_as_advanced(
CPLEX_LIBRARIES
CPLEX_INCLUDE_DIR)

1
googletest

@ -0,0 +1 @@
Subproject commit aa148eb2b7f70ede0eb10de34b6254826bfb34f4

2
infinity/benchmark/CMakeLists.txt
Unescape View File

@ -0,0 +1,2 @@
add_executable(infinity-benchmark.run src/main.c)
target_link_libraries (infinity-benchmark.run LINK_PUBLIC lifting_static infinity_static multirow_static)

467
infinity/benchmark/src/main.c
Unescape View File

@ -0,0 +1,467 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define _XOPEN_SOURCE 500
#include <getopt.h>
#include <unistd.h>
#include <time.h>
#include <multirow/cg.h>
#include <multirow/mir.h>
#include <multirow/stats.h>
#include <multirow/util.h>
#include <infinity/greedy.h>
int ENABLE_LIFTING = 0;
int MIN_N_ROWS = 2;
int MAX_N_ROWS = 2;
int DUMP_CUT = 0;
int DUMP_CUT_N = 0;
int GENERATE_MIR = 0;
int GENERATE_GREEDY = 0;
int KEEP_INTEGRALITY = 0;
int N_ROUNDS = 1;
char BASIS_FILENAME[1000] = {0};
char PROBLEM_FILENAME[1000] = {0};
char KNOWN_SOLUTION_FILENAME[1000] = {0};
char OUTPUT_SOLUTION_FILENAME[1000] = {0};
char OUTPUT_BASIS_FILENAME[1000] = {0};
char LOG_FILENAME[1000] = {0};
char STATS_FILENAME[1000] = {0};
FILE *LOG_FILE;
int BOOST_VAR = -1;
double BOOST_FACTOR = 0.01;
#define OPTION_WRITE_BASIS 1000
#define OPTION_WRITE_SOLUTION 1001
static const struct option options_tab[] =
{
{"help", no_argument, 0, 'h'},
{"problem", required_argument, 0, 'p'},
{"solution", required_argument, 0, 'x'},
{"mir", no_argument, 0, 'm'},
{"greedy", no_argument, 0, 'g'},
{"rounds", required_argument, 0, 'r'},
{"keep-integrality", no_argument, 0, 'k'},
{"write-solution", required_argument, 0, OPTION_WRITE_SOLUTION},
{"write-basis", required_argument, 0, OPTION_WRITE_BASIS},
{"basis", required_argument, 0, 'b'},
{"log", required_argument, 0, 'l'},
{"stats", required_argument, 0, 's'},
{"boost", required_argument, 0, 't'},
{"lift", no_argument, 0, 'i'},
{"rows", required_argument, 0, 'a'},
{0, 0, 0, 0}
};
static void print_usage(char **argv)
{
printf("Usage: %s [OPTION]...\n", argv[0]);
printf("Solves the given MILP using specified cutting-planes.\n\n");
printf("Parameters:\n");
printf("%4s %-20s %s\n", "-b", "--basis=FILE",
"BAS file containing an optimal basis for the linear relaxation of "
"the problem");
printf("%4s %-20s %s\n", "-g", "--greedy",
"generate greedy intersection cuts");
printf("%4s %-20s %s\n", "-k", "--keep-integrality",
"do not relax integrality of variables");
printf("%4s %-20s %s\n", "-l", "--log=FILE",
"write log to the specified file");
printf("%4s %-20s %s\n", "-m", "--mir", "generate MIR cuts");
printf("%4s %-20s %s\n", "-p", "--problem=FILE", "problem to be solved");
printf("%4s %-20s %s\n", "-r", "--rounds=NUM",
"number of rounds of cutting planes to add");
printf("%4s %-20s %s\n", "-s", "--stats=FILE",
"write statistics to the specified file");
printf("%4s %-20s %s\n", "-x", "--solution=FILE",
"known integral solution (to check the validity of the cuts)");
printf("%4s %-20s %s\n", "", "--write-solution=FILE",
"write solution found at the end of the procedure to given file");
printf("%4s %-20s %s\n", "", "--write-basis=FILE",
"write optimal LP basis to given file");
printf("%4s %-20s %s\n", "", "--lift", "enable trivial lifting");
printf("%4s %-20s %s\n", "", "--rows=N",
"generate multi-row cuts from up to N rows");
}
static int parse_args(int argc,
char **argv)
{
int rval = 0;
opterr = 0;
while (1)
{
int c = 0;
int option_index = 0;
c = getopt_long(argc, argv, "l:p:s:x:hmgr:kb:t:ia:", options_tab,
&option_index);
if (c < 0) break;
switch (c)
{
case 'l':
strcpy(LOG_FILENAME, optarg);
break;
case 'b':
strcpy(BASIS_FILENAME, optarg);
break;
case 'g':
GENERATE_GREEDY = 1;
break;
case 'k':
KEEP_INTEGRALITY = 1;
break;
case 'm':
GENERATE_MIR = 1;
break;
case 't':
BOOST_VAR = atoi(optarg);
break;
case 'r':
N_ROUNDS = atoi(optarg);
break;
case 'a':
MAX_N_ROWS = atoi(optarg);
break;
case OPTION_WRITE_SOLUTION:
strcpy(OUTPUT_SOLUTION_FILENAME, optarg);
break;
case OPTION_WRITE_BASIS:
strcpy(OUTPUT_BASIS_FILENAME, optarg);
break;
case 'p':
strcpy(PROBLEM_FILENAME, optarg);
break;
case 's':
strcpy(STATS_FILENAME, optarg);
break;
case 'x':
strcpy(KNOWN_SOLUTION_FILENAME, optarg);
break;
case 'h':
print_usage(argv);
exit(0);
case 'i':
ENABLE_LIFTING = 1;
break;
case ':':
fprintf(stderr, "%s: option '-%c' requires an argument\n",
argv[0], optopt);
rval = 1;
goto CLEANUP;
case '?':
default:
fprintf(stderr, "%s: option '-%c' is invalid\n", argv[0],
optopt);
rval = 1;
goto CLEANUP;
}
}
if ((strlen(PROBLEM_FILENAME) == 0))
{
fprintf(stderr, "You must specify the problem.\n");
rval = 1;
}
if (KEEP_INTEGRALITY && (GENERATE_GREEDY + GENERATE_MIR > 0))
{
fprintf(stderr, "Cutting planes cannot be added when integrality is "
"kept\n");
rval = 1;
}
if (N_ROUNDS < 1)
{
fprintf(stderr, "Invalid number of rounds.\n");
rval = 1;
}
CLEANUP:
if (rval)
fprintf(stderr, "Try '%s --help' for more information.\n", argv[0]);
return rval;
}
void print_header(int argc,
char *const *argv)
{
char hostname[3000];
gethostname(hostname, 1024);
time_t now;
time(&now);
struct tm *ttm = localtime(&now);
time_printf("multirow\n");
time_printf("%s %04d-%02d-%02d %02d:%02d\n", hostname, ttm->tm_year + 1900,
ttm->tm_mon + 1, ttm->tm_mday, ttm->tm_hour, ttm->tm_min);
time_printf("Compile-time parameters:\n");
time_printf(" EPSILON: %e\n", EPSILON);
time_printf(" GREEDY_BIG_E: %e\n", GREEDY_BIG_E);
time_printf(" GREEDY_MAX_GAP: %e\n", GREEDY_MAX_GAP);
char cmdline[5000] = {0};
for (int i = 0; i < argc; i++)
sprintf(cmdline + strlen(cmdline), "%s ", argv[i]);
time_printf("Command line arguments:\n");
time_printf(" %s\n", cmdline);
}
int main(int argc,
char **argv)
{
int rval = 0;
double *x = 0;
struct CG *cg = 0;
struct LP lp;
char *column_types = 0;
rval = parse_args(argc, argv);
if (rval) return 1;
if (LOG_FILENAME[0])
{
LOG_FILE = fopen(LOG_FILENAME, "w");
abort_if(!LOG_FILE, "could not open log file");
}
if_info_level
{
print_header(argc, argv);
}
STATS_init();
STATS_set_input_filename(PROBLEM_FILENAME);
progress_title(PROBLEM_FILENAME);
rval = LP_open(&lp);
abort_if(rval, "LP_open failed");
rval = LP_create(&lp, "multirow");
abort_if(rval, "LP_create failed");
rval = LP_read_problem(&lp, PROBLEM_FILENAME);
abort_if(rval, "LP_read_problem failed");
int ncols = LP_get_num_cols(&lp);
int nrows = LP_get_num_rows(&lp);
log_info(" %d rows, %d cols\n", nrows, ncols);
x = (double *) malloc(ncols * sizeof(double));
abort_if(!x, "could not allocate x");
if (!KEEP_INTEGRALITY)
{
column_types = (char *) malloc(ncols * sizeof(char));
abort_if(!column_types, "could not allocate column_types");
log_info("Storing column types...\n");
rval = LP_get_column_types(&lp, column_types);
abort_if(rval, "LP_get_column_types failed");
log_info("Relaxing integrality...\n");
rval = LP_relax(&lp);
abort_if(rval, "LP_relax failed");
log_info("Disabling presolve...\n");
LP_disable_presolve(&lp);
}
if(BASIS_FILENAME[0])
{
rval = LP_read_basis(&lp, BASIS_FILENAME);
abort_if(rval, "LP_read_basis failed");
}
log_info("Optimizing...\n");
int infeasible;
rval = LP_optimize(&lp, &infeasible);
abort_if(rval, "LP_optimize failed");
double cost = 0, xboost = 0;
if(BOOST_VAR > 0)
{
rval = LP_get_x(&lp, x);
abort_if(rval, "LP_get_x failed");
xboost = x[BOOST_VAR];
for (int i = 0; i < ncols; i++)
log_info("x[%3d] = %12.8lf\n", i, x[i]);
}
double obj;
rval = LP_get_obj_val(&lp, &obj);
abort_if(rval, "LP_get_obj_val failed");
log_info(" opt = %lf\n", obj);
STATS_set_obj_value(obj);
STATS_finish_round();
if(BOOST_VAR >= 0)
log_info("Boosting variable %d by %.2lf\n", BOOST_VAR, BOOST_FACTOR);
if(OUTPUT_BASIS_FILENAME[0])
{
rval = LP_write_basis(&lp, OUTPUT_BASIS_FILENAME);
abort_if(rval, "LP_write_basis failed");
}
if(GENERATE_MIR || GENERATE_GREEDY)
{
cg = (struct CG *) malloc(sizeof(struct CG));
abort_if(!cg, "could not allocate cg");
log_info("Reading tableau rows...\n");
rval = CG_init(&lp, column_types, cg);
abort_if(rval, "CG_init failed");
if (strlen(KNOWN_SOLUTION_FILENAME) > 0)
{
if(access(KNOWN_SOLUTION_FILENAME, F_OK) != -1)
{
rval = LP_read_solution(&lp, KNOWN_SOLUTION_FILENAME, x);
abort_if(rval, "LP_read_solution failed");
rval = CG_set_integral_solution(cg, x);
abort_if(rval, "CG_set_integral_solution failed");
}
else
{
log_error("ERROR: Could not read solution file!\n");
}
}
rval = LP_get_x(&lp, x);
abort_if(rval, "LP_get_x failed");
rval = CG_set_basic_solution(cg, x);
abort_if(rval, "CG_set_basic_solution failed");
if (GENERATE_MIR)
{
log_info("Adding MIR cuts...\n");
rval = CG_add_single_row_cuts(
cg,
(SingleRowGeneratorCallback)
MIR_generate_cut);
abort_if(rval, "CG_add_single_row_cuts failed");
log_info("Optimizing...\n");
rval = LP_optimize(&lp, &infeasible);
abort_if(rval, "LP_optimize failed");
rval = LP_get_obj_val(&lp, &obj);
abort_if(rval, "LP_get_obj_val failed");
log_info(" opt = %lf\n", obj);
STATS_set_obj_value(obj);
STATS_finish_round();
}
if (GENERATE_GREEDY)
{
for(int k = MIN_N_ROWS; k <= MAX_N_ROWS; k++)
{
log_info("Adding greedy intersection cuts (%d rows)...\n", k);
rval = CG_add_multirow_cuts(cg, k, (MultirowGeneratorCallback)
GREEDY_generate_cut);
abort_if(rval, "CG_add_multirow_cuts failed");
}
log_info("Optimizing...\n");
rval = LP_optimize(&lp, &infeasible);
abort_if(rval, "LP_optimize failed");
rval = LP_get_obj_val(&lp, &obj);
abort_if(rval, "LP_get_obj_val failed");
log_info(" opt = %lf\n", obj);
STATS_set_obj_value(obj);
STATS_finish_round();
}
CG_free(cg);
}
if(BOOST_VAR > 0)
{
FILE *boost = fopen("boost.csv", "a");
fprintf(boost, "%d, %.8lf, %.8lf, %.8lf\n", BOOST_VAR, obj, cost,
xboost);
}
if(OUTPUT_SOLUTION_FILENAME[0])
{
rval = LP_write_solution(&lp, OUTPUT_SOLUTION_FILENAME);
abort_if(rval, "LP_write_solution failed");
}
if(STATS_FILENAME[0])
{
log_info("Writing stats to file %s...\n", STATS_FILENAME);
rval = STATS_print_yaml(STATS_FILENAME);
abort_if(rval, "STATS_print_yaml failed");
}
CLEANUP:
if (LOG_FILE) fclose(LOG_FILE);
if (x) free(x);
if (column_types) free(column_types);
LP_free(&lp);
return rval;
}

20
infinity/library/CMakeLists.txt
Unescape View File

@ -0,0 +1,20 @@
set(COMMON_SOURCES
src/greedy-2d.c
src/greedy-nd.c
src/greedy-bsearch.c
src/greedy.c
include/infinity/greedy-2d.h
include/infinity/greedy-nd.h
include/infinity/greedy-bsearch.h
include/infinity/greedy.h)
set(TEST_SOURCES
tests/greedy-2d-test.cpp
tests/greedy-nd-test.cpp)
add_library(infinity_static ${COMMON_SOURCES})
set_target_properties(infinity_static PROPERTIES OUTPUT_NAME lifting)
target_include_directories (infinity_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
add_executable(infinity-test.run ${COMMON_SOURCES} ${TEST_SOURCES})
target_link_libraries(infinity-test.run gtest_main multirow_static lifting_static)

37
infinity/library/include/infinity/greedy-2d.h
Unescape View File

@ -0,0 +1,37 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_GREEDY_2D_H
#define MULTIROW_GREEDY_2D_H
int GREEDY_2D_bound(const double *rays,
const double *bounds,
int nrays,
const double *f,
const double *p,
double *epsilon,
double *v1,
double *v2,
int *index1,
int *index2);
int GREEDY_2D_generate_cut(const double *rays,
int nrays,
const double *f,
double *bounds);
#endif //MULTIROW_GREEDY_2D_H

33
infinity/library/include/infinity/greedy-bsearch.h
Unescape View File

@ -0,0 +1,33 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_GREEDY_BSEARCH_H
#define MULTIROW_GREEDY_BSEARCH_H
int create_sfree_mip(int nrows,
int nrays,
const double *f,
const double *rays,
const double *bounds,
double e,
struct LP *lp);
int GREEDY_BSEARCH_compute_bounds(int nrows,
int nrays,
const double *f,
const double *rays,
double *bounds);
#endif //MULTIROW_GREEDY_BSEARCH_H

105
infinity/library/include/infinity/greedy-nd.h
Unescape View File

@ -0,0 +1,105 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#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,
const double *rays,
const double *beta,
struct LP *lp);
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 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 GREEDY_ND_generate_cut(int nrows,
int nrays,
const double *f,
const double *rays,
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 GREEDY_ND_cone_bound(int nrows,
int nrays,
const double *f,
const double *rays,
const int *rx,
const double *x,
const double *beta,
double *epsilon);
int GREEDY_ND_find_violated_cone(int nrows,
int nrays,
const double *f,
const double *rays,
const double *x,
const double *beta,
double epsilon,
int *rx,
double *sbar,
int *violated_found);
int GREEDY_ND_find_tight_rays(int nrows,
int nrays,
const double *f,
const double *rays,
const double *x,
const double *beta,
double epsilon,
int *tx);
int GREEDY_ND_scale_to_ahull(int nrows,
int nrays,
const double *rays,
const int *rx,
const double *beta,
double epsilon,
const double *d,
double *alpha);
#endif //MULTIROW_GREEDY_ND_H

31
infinity/library/include/infinity/greedy.h
Unescape View File

@ -0,0 +1,31 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_GREEDY_H
#define MULTIROW_GREEDY_H
int GREEDY_write_sage_file(int nrows,
int nrays,
const double *f,
const double *rays,
const double *bounds,
const char *filename);
int GREEDY_generate_cut(int nrows,
struct Row **rows,
const char *column_types,
struct Row *cut);
#endif //MULTIROW_GREEDY_H

738
infinity/library/src/greedy-2d.c
Unescape View File

@ -0,0 +1,738 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define _GNU_SOURCE
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <multirow/geometry.h>
#include <multirow/double.h>
#include <multirow/util.h>
#include <multirow/rational.h>
#include <infinity/greedy-2d.h>
static int get_bounding_box(int nrows,
int nrays,
const double *rays,
const double *bounds,
double epsilon,
double *lb,
double *ub)
{
for(int i = 0; i < nrows; i++)
{
ub[i] = 0;
lb[i] = 0;
}
for (int i=0; i < nrays; i++)
{
double e = fmin(epsilon, bounds[i]);
for(int j = 0; j < nrows; j++)
{
double rij = rays[nrows * i + j] * e;
lb[j] = fmin(lb[j], floor(rij));
ub[j] = fmax(ub[j], ceil(rij));
}
}
return 0;
}
struct LatticeSequence
{
int square;
int direction;
int steps;
int eol;
int i;
int j;
};
static void lattice_sequence_init(struct LatticeSequence *seq)
{
seq->steps = 0;
seq->square = 0;
seq->direction = 0;
seq->i = 0;
seq->j = 0;
seq->eol = 0;
}
static int next_lattice_point(struct LatticeSequence *seq,
const double *lb,
const double *ub)
{
int rval = 0;
seq->eol = 1;
for(int k = 0; k < 8 * (seq->square+1); k++)
{
if(seq->steps > 0)
{
seq->steps--;
switch(seq->direction)
{
case 3:
seq->i++;
break;
case 2:
seq->j--;
break;
case 1:
seq->i--;
break;
case 0:
seq->j++;
break;
}
}
if(seq->steps == 0)
{
if(seq->direction > 0)
{
seq->direction--;
seq->steps = 2 * seq->square;
}
else
{
seq->square++;
seq->direction = 3;
seq->steps = 2 * seq->square;
seq->i--;
seq->j++;
}
}
if(seq->i >= lb[0] && seq->i <= ub[0] && seq->j >= lb[1] &&
seq->j <= ub[1])
{
seq->eol = 0;
break;
}
}
CLEANUP:
return rval;
}
/*
* Returns lambda1 and lambda2 such that p = r1 * lambda1 + r2 * lambda2
*/
static void find_linear_combination(const double *r1,
const double *r2,
const double *p,
double *lambda1,
double *lambda2)
{
int rval = 0;
double den = (r1[0] * r2[1] - r2[0] * r1[1]);
if(DOUBLE_iszero(den)) den = 0.0;
*lambda1 = (r2[1] * p[0] - r2[0] * p[1]) / den;
*lambda2 = (-r1[1] * p[0] + r1[0] * p[1]) / den;
}
static int generate_split(const double *f,
const double *d,
double *pi,
double *pi_zero,
long max_den)
{
int rval = 0;
Rational d1, d2;
double m;
rval = DOUBLE_to_rational(d[0], 10, d1);
abort_if(rval, "DOUBLE_to_rational failed");
rval = DOUBLE_to_rational(d[1], 10, d2);
abort_if(rval, "DOUBLE_to_rational failed");
m = lcm(d1->den, d2->den);
d1->den = m / d1->den;
d1->num *= d1->den;
d2->den = m / d2->den;
d2->num *= d2->den;
m = gcd(d1->num, d2->num);
if(m != 0)
{
d1->num /= m;
d2->num /= m;
pi[0] = (double) d2->num;
pi[1] = - ((double) d1->num);
*pi_zero = floor(pi[0] * f[0] + pi[1] * f[1]);
}
else
{
pi[0] = pi[1] = INFINITY;
*pi_zero = INFINITY;
}
CLEANUP:
return rval;
}
/*
* Receives a list of rays r1,...,rn and a point p. Returns i1,i2 such that
* p belongs to cone(r_i1,r_i2). Also returns lambda1, lambda2 such that
* p = r_i1 * lambda1 + r_i2 * lambda2.
*
* The rays must be sorted in clockwise order.
*/
static int find_containing_cone(const double *rays,
const int nrays,
const double *p,
int *index1,
int *index2,
double *lambda1,
double *lambda2)
{
int rval = 0;
int i1, i2;
for (i1 = 0; i1 < nrays; i1++)
{
i2 = (i1 + 1) % nrays;
const double *r1 = &rays[i1 * 2];
const double *r2 = &rays[i2 * 2];
double at1 = atan2(r1[1], r1[0]);
double at2 = atan2(r2[1], r2[0]) - at1;
if (at2 > 0) at2 -= 2 * M_PI;
if (at2 <= - M_PI)
{
log_verbose(" discarding obtuse rays\n");
log_verbose(" r1=%.12lf %.12lf\n", r1[0], r1[1]);
log_verbose(" r2=%.12lf %.12lf\n", r2[0], r2[1]);
continue;
}
find_linear_combination(r1, r2, p, lambda1, lambda2);
log_verbose(" r1=%.12lf %.12lf\n", r1[0], r1[1]);
log_verbose(" r2=%.12lf %.12lf\n", r2[0], r2[1]);
log_verbose(" %.8lf %.8lf\n", *lambda1, *lambda2);
if(DOUBLE_iszero(*lambda1)) *lambda1 = 0.0;
if(DOUBLE_iszero(*lambda2)) *lambda2 = 0.0;
if ((*lambda1) >= 0 && (*lambda2) >= 0 && (*lambda1) <= 1e9 &&
(*lambda2) <= 1e9)
{
log_verbose(" match!\n");
break;
}
}
if (i1 == nrays)
i1 = i2 = -1;
*index1 = i1;
*index2 = i2;
CLEANUP:
return rval;
}
/*
* Find lambda such that p lies on the line connecting lambda * r1 and lambda * r2
*/
static int scale_cone_to_point(const double *r1,
const double *r2,
const double *p,
double *lambda)
{
int rval = 0;
double a = r1[0], b = r1[1];
double c = r2[0], d = r2[1];
double den = (b * c - a * d);
//abort_iff(fabs(den) < 1e-9, "rays cannot be parallel (den=%.12lf)", den);
*lambda = p[0] * (b - d) - p[1] * (a - c);
*lambda /= den;
CLEANUP:
return rval;
}
/*
* Find lambda such that p lies in the line connecting r1 and lambda * r2
*/
static int shear_cone_to_point(const double *r1,
const double *r2,
const double *p,
double *lambda)
{
int rval = 0;
double a = r1[0], b = r1[1];
double c = r2[0], d = r2[1];
double den = d * (p[0] - a) - c * (p[1] - b);
*lambda = b * p[0] - a * p[1];
*lambda /= den;
CLEANUP:
return rval;
}
static int scale_to_chull(double *rays, int nrays, double *scale)
{
int rval = 0;
double *rays_extended = 0;
double *vertices = 0;
int nvertices;
rays_extended = (double*) malloc(2 * (nrays + 1) * sizeof(double));
vertices = (double*) malloc(2 * (nrays + 1) * sizeof(double));
abort_if(!rays_extended, "could not allocate rays_extended");
abort_if(!vertices, "could not allocate vertices");
memcpy(rays_extended, rays, 2 * nrays * sizeof(double));
rays_extended[2*nrays] = 0;
rays_extended[2*nrays + 1] = 0;
rval = chull_2d(rays_extended, nrays + 1, vertices, &nvertices);
abort_if(rval, "chull_2d failed");
for(int i = 0; i < nrays; i++)
scale[i] = 1.0;
log_verbose(" convex hull:\n");
for(int i = 0; i < nvertices; i++)
{
log_verbose(" v%-3d: %20.8lf %20.8lf\n", i, vertices[2 * i],
vertices[2 * i + 1]);
}
log_verbose(" rays:\n");
for(int i = 0; nrays >= 3 && i < nrays; i++)
{
int i1, i2;
double lambda1, lambda2, mu;
rval = find_containing_cone(vertices, nvertices, &rays[2*i], &i1, &i2,
&lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
log_verbose("%.8lf %.8lf\n", lambda1, lambda2);
if(i1 < 0 || i2 < 0) continue;
if(DOUBLE_iszero(lambda1))
{
mu = lambda2;
}
else if(DOUBLE_iszero(lambda2))
{
mu = lambda1;
}
else
{
rval = scale_vector_to_line(&vertices[2*i1], &vertices[2*i2],
&rays[2*i], &mu);
abort_if(rval, "scale_vector_to_line failed");
}
abort_if(!isfinite(mu), "mu should be finite");
//log_verbose(" r%-3d: %.2lf %.2lf %.2lf\n", i, rays[2*i], rays[2*i+1], mu);
rays[2*i] *= mu;
rays[2*i+1] *= mu;
scale[i] = mu;
log_verbose(" r%-3d: %20.12lf %20.12lf (scale %.8lf)\n", i, rays[2*i],
rays[2*i+1], scale[i]);
}
CLEANUP:
if(rays_extended) free(rays_extended);
if(vertices) free(vertices);
return rval;
}
#ifndef TEST_SOURCE
int GREEDY_2D_bound(const double *rays,
const double *bounds,
int nrays,
const double *f,
const double *p,
double *epsilon,
double *v1,
double *v2,
int *index1,
int *index2)
{
int rval = 0;
int i1, i2, iexact = -1;
double e1, e2;
const double *r1, *r2;
double lambda1, lambda2;
double pp[2] = { p[0] - f[0], p[1] - f[1] };
rval = find_containing_cone(rays, nrays, pp, &i1, &i2, &lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
if(i1 < 0 || i2 < 0)
{
log_verbose(" no cone\n");
*epsilon = INFINITY;
goto CLEANUP;
}
if(DOUBLE_iszero(lambda1)) iexact = i2;
if(DOUBLE_iszero(lambda2)) iexact = i1;
if(iexact >= 0)
{
int inext = (iexact + 1) % nrays;
int iprev = (iexact + (nrays-1)) % nrays;
double mu1, mu2;
find_linear_combination(&rays[2 * iprev], &rays[2 * inext],
&rays[2 * iexact], &mu1, &mu2);
log_verbose(" mu=%.12lf %.12lf\n", mu1, mu2);
int should_enlarge_cone = 1;
if(!DOUBLE_leq(bounds[iexact], lambda1+lambda2)) should_enlarge_cone = 0;
if(!DOUBLE_geq(mu1, 0) || !DOUBLE_geq(mu2, 0)) should_enlarge_cone = 0;
if(isinf(mu1) || isinf(mu2)) should_enlarge_cone = 0;
if(fabs(mu1) > 1e9 || fabs(mu2) > 1e9) should_enlarge_cone = 0;
if(should_enlarge_cone)
{
i1 = iprev;
i2 = inext;
}
else
{
v1[0] = v1[1] = INFINITY;
v2[0] = v2[1] = INFINITY;
*epsilon = (lambda1 + lambda2);
if(DOUBLE_leq(bounds[iexact], *epsilon))
*epsilon = INFINITY;
*index1 = *index2 = iexact;
goto CLEANUP;
}
}
if(bounds[i1] > bounds[i2])
{
swap(i1, i2, int);
swap(lambda1, lambda2, double);
}
*index1 = i1;
*index2 = i2;
r1 = &rays[2 * i1];
r2 = &rays[2 * i2];
log_verbose(" r%-4d %20.12lf %20.12lf\n", i1, r1[0], r1[1]);
log_verbose(" r%-4d %20.12lf %20.12lf\n", i2, r2[0], r2[1]);
log_verbose(" pp %20.12lf %20.12lf\n", pp[0], pp[1]);
log_verbose(" lambda %20.12lf %20.12lf\n", lambda1, lambda2);
double r1bound[2] = { r1[0] * bounds[i1], r1[1] * bounds[i1] };
rval = scale_cone_to_point(r1, r2, pp, &e1);
abort_if(rval, "scale_cone_to_point failed");
rval = shear_cone_to_point(r1bound, r2, pp, &e2);
abort_if(rval, "scale_cone_to_point failed");
log_verbose(" e1=%20.12lf\n", e1);
log_verbose(" e2=%20.12lf\n", e2);
log_verbose(" b1=%20.12lf\n", bounds[i1]);
log_verbose(" b2=%20.12lf\n", bounds[i2]);
switch(DOUBLE_cmp(e1, bounds[i1]))
{
case -1:
*epsilon = e1;
v1[0] = r1[0] * e1;
v1[1] = r1[1] * e1;
v2[0] = r2[0] * e1;
v2[1] = r2[1] * e1;
break;
case 0:
case 1:
if(DOUBLE_geq(e2, bounds[i2]))
{
*epsilon = INFINITY;
}
else
{
*epsilon = e2;
v1[0] = r1[0] * bounds[i1];
v1[1] = r1[1] * bounds[i1];
v2[0] = r2[0] * e2;
v2[1] = r2[1] * e2;
}
}
CLEANUP:
return rval;
}
int GREEDY_2D_generate_cut(const double *original_rays,
const int nrays,
const double *f,
double *bounds)
{
log_verbose("GREEDY_2D_generate_cut\n");
int rval = 0;
int count = 0;
double *scale = 0;
double *rays = 0;
double lb[2], ub[2];
for (int i = 0; i < nrays; i++)
bounds[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, original_rays, 2 * nrays * sizeof(double));
rval = scale_to_chull(rays, nrays, scale);
abort_if(rval, "scale_to_chull failed");
long seq_count = 0;
while(1)
{
log_verbose(" starting iteration %d...\n", count);
abort_if(count++ > 2 * nrays, "infinite loop");
rval = get_bounding_box(2, nrays, rays, bounds, 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]);
int best_i1, best_i2;
double best_v1[2];
double best_v2[2];
double best_p[2];
double best_epsilon = INFINITY;
struct LatticeSequence seq;
lattice_sequence_init(&seq);
while(!seq.eol)
{
seq_count++;
double p[2] = { seq.i, seq.j };
double v1[2], v2[2], epsilon;
int i1, i2;
log_verbose(" p=%.2lf %.2lf\n", p[0], p[1]);
rval = GREEDY_2D_bound(rays, bounds, nrays, f, p, &epsilon, v1, v2,
&i1, &i2);
abort_if(rval, "GREEDY_2D_bound failed");
log_verbose(" epsilon=%.2lf\n", epsilon);
if(epsilon >= 0 && epsilon < best_epsilon)
{
log_verbose(" found smaller epsilon: %.8lf\n", epsilon);
rval = get_bounding_box(2, nrays, rays, bounds, epsilon, lb, ub);
abort_if(rval, "get_bounding_box failed");
log_verbose(" p=%.2lf %.2lf\n", p[0], p[1]);
log_verbose(" box=[%.2lf %.2lf] [%.2lf %.2lf]\n", lb[0], ub[0],
lb[1], ub[1]);
log_verbose(" v1=%12.8lf %12.8lf\n", v1[0], v1[1]);
log_verbose(" v2=%12.8lf %12.8lf\n", v2[0], v2[1]);
log_verbose(" rays=[%d %d]\n", i1, i2);
best_epsilon = epsilon;
best_v1[0] = v1[0];
best_v1[1] = v1[1];
best_v2[0] = v2[0];
best_v2[1] = v2[1];
best_p[0] = p[0];
best_p[1] = p[1];
best_i1 = i1;
best_i2 = i2;
}
next_lattice_point(&seq, lb, ub);
if(seq_count > MAX_LATTICE_POINTS)
{
rval = ERR_MIP_TIMEOUT;
goto CLEANUP;
}
}
if(isinf(best_epsilon))
{
log_verbose(" best_epsilon is infinity\n");
break;
}
log_verbose(" updating bounds\n");
if(isinf(best_v1[0]))
{
bounds[best_i1] = best_epsilon;
log_verbose(" bounds[%d]=%.8lf (exact)\n", best_i1, best_epsilon);
}
else
{
log_verbose(" v1=%.8lf %.8lf\n", best_v1[0], best_v1[1]);
log_verbose(" v2=%.8lf %.8lf\n", best_v2[0], best_v2[1]);
log_verbose(" i=%d %d\n", best_i1, best_i2);
for (int i = 0; i < nrays; i++)
{
double lambda;
rval = scale_vector_to_line(best_v1, best_v2, &rays[2 * i], &lambda);
abort_if(rval, "scale_vector_to_line failed");
if(!DOUBLE_geq(lambda, 0)) continue;
bounds[i] = fmin(bounds[i], lambda);
log_verbose(" bounds[%d]=%.8lf\n", i, bounds[i]);
}
}
//if(count > 0)
//{
// for (int i = 0; i < nrays; i++)
// bounds[i] = fmin(bounds[i], best_epsilon);
// break;
//}
int is_split;
for (int k = 0; k < nrays; k++)
{
if(bounds[k] < 100) continue;
is_split = 1;
double *split_direction = &rays[2 * k];
log_verbose(" split_direction=%.2lf %.2lf\n", split_direction[0],
split_direction[1]);
double pi[2];
double pi_zero;
rval = generate_split(f, split_direction, pi, &pi_zero, 10);
abort_if(rval, "generate_split failed");
log_verbose(" pi=%.2lf %.2lf\n", pi[0], pi[1]);
log_verbose(" pi_zero=%.2lf\n", pi_zero);
if(isinf(pi[0]))
{
is_split = 0;
break;
}
double lhs;
// reject splits that have f on the boundary
lhs = f[0] * pi[0] + f[1] * pi[1];
if(DOUBLE_eq(pi_zero, lhs) || DOUBLE_eq(lhs, pi_zero+1))
{
log_verbose(" split rejected\n");
is_split = 0;
}
for (int i = 0; i < nrays && is_split; i++)
{
const double *r = &rays[2 * i];
lhs = (f[0] + r[0] * bounds[i]) * pi[0];
lhs += (f[1] + r[1] * bounds[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]);
is_split = 0;
}
}
if(is_split)
{
log_verbose(" split confirmed. stopping.\n");
break;
}
}
if(is_split) break;
}
for(int i=0; i<nrays; i++)
bounds[i] *= scale[i];
CLEANUP:
if(scale) free(scale);
if(rays) free(rays);
return rval;
}
#endif // TEST_SOURCE

246
infinity/library/src/greedy-bsearch.c
Unescape View File

@ -0,0 +1,246 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <stdlib.h>
#include <multirow/cg.h>
#include <multirow/double.h>
#include <multirow/geometry.h>
#include <multirow/lp.h>
#include <multirow/util.h>
#include <infinity/greedy-bsearch.h>
int create_sfree_mip(int nrows,
int nrays,
const double *f,
const double *rays,
const double *bounds,
double e,
struct LP *lp)
{
int rval = 0;
double rhs;
char sense;
int rmatbeg = 0;
int* rmatind = 0;
double *rmatval = 0;
rmatind = (int *) malloc((nrows + nrays) * sizeof(int));
rmatval = (double *) malloc((nrows + nrays) * sizeof(double));
abort_if(!rmatind, "could not allocate rmatind");
abort_if(!rmatval, "could not allocate rmatval");
rval = LP_create(lp, "greedy");
abort_if(rval, "LP_create failed");
// create x (basic) variables
for (int i = 0; i < nrows; i++)
{
rval = LP_new_col(lp, 0, -MILP_INFINITY, MILP_INFINITY, 'I');
abort_if(rval, "LP_new_col failed");
}
// create s (non-basic) variables
for (int i = 0; i < nrays; i++)
{
rval = LP_new_col(lp, 1.0, 0, MILP_INFINITY, 'C');
abort_if(rval, "LP_new_col failed");
}
// add constraint \sum_{i=1}^m s_i \leq 1
sense = 'L';
rhs = 1.0;
for (int i = 0; i < nrays; i++)
{
rmatind[i] = i + nrows;
rmatval[i] = 1.0;
}
rval = LP_add_rows(lp, 1, nrays, &rhs, &sense, &rmatbeg, rmatind, rmatval);
abort_if(rval, "LP_add_rows failed");
// add constraints x_i - \sum_{j=1}^m min{e,e_j} s_j R_ji = f_i
for (int i = 0; i < nrows; i++)
{
int k = 0;
sense = 'E';
rhs = f[i];
rmatind[k] = i;
rmatval[k] = 1.0;
k++;
for (int j = 0; j < nrays; j++)
{
rmatind[k] = j + nrows;
rmatval[k] = -rays[nrows * j + i] * fmin(e, bounds[j]);
k++;
}
rval = LP_add_rows(lp, 1, nrays + 1, &rhs, &sense, &rmatbeg, rmatind,
rmatval);
abort_if(rval, "LP_add_rows failed");
}
CLEANUP:
if (rmatind) free(rmatind);
if (rmatval) free(rmatval);
return rval;
}
int GREEDY_BSEARCH_compute_bounds(int nrows,
int nrays,
const double *f,
const double *rays,
double *bounds)
{
int rval = 0;
struct LP lp;
double e_upper = 2 * GREEDY_BIG_E;
double e_lower = 0.0;
int cplex_count = 0;
double cplex_time = 0;
int iteration_count = 0;
double *x = 0;
x = (double *) malloc((nrays + nrows) * sizeof(double));
abort_if(!x, "could not allocate x");
for (int i = 0; i < nrays; i++)
bounds[i] = GREEDY_BIG_E;
for (int it = 0;; it++)
{
abort_if(it > 2*nrays, "stuck in an infinite loop");
log_verbose("Starting iteration %d...\n", it);
iteration_count++;
int solution_found = 0;
int inner_count = 0;
while (fabs(e_upper - e_lower) > GREEDY_MAX_GAP)
{
inner_count++;
double e = (e_upper + e_lower) / 2;
log_verbose(" e=%.12lf\n", e);
rval = LP_open(&lp);
abort_if(rval, "LP_open failed");
rval = create_sfree_mip(nrows, nrays, f, rays, bounds, e, &lp);
abort_if(rval, "create_sfree_mip failed");
if_verbose_level
{
rval = LP_write(&lp, "greedy.lp");
abort_if(rval, "LP_write failed");
}
int infeasible;
cplex_count++;
double initial_time = get_user_time();
log_verbose(" Optimizing...\n");
rval = LP_optimize(&lp, &infeasible);
if (rval)
{
// Workaround for CPLEX bug. If CPLEX tell us that this problem
// is unbounded, we disable presolve and try again.
LP_free(&lp);
LP_open(&lp);
rval = create_sfree_mip(nrows, nrays, f, rays, bounds, e, &lp);
abort_if(rval, "create_sfree_mip failed");
LP_disable_presolve(&lp);
rval = LP_optimize(&lp, &infeasible);
abort_if(rval, "LP_optimize failed");
}
cplex_time += get_user_time() - initial_time;
if (infeasible)
{
e_lower = e;
log_verbose(" infeasible\n");
if (e > GREEDY_BIG_E-1)
{
LP_free(&lp);
goto OUT;
}
}
else
{
log_verbose(" feasible\n");
e_upper = e;
solution_found = 1;
rval = LP_get_x(&lp, x);
abort_if(rval, "LP_get_x failed");
}
LP_free(&lp);
}
if (solution_found)
{
for (int j = 0; j < nrays; j++)
{
if (!DOUBLE_geq(x[nrows + j], 0.001)) continue;
bounds[j] = fmin(bounds[j] * 0.99, e_lower * 0.99);
}
}
log_verbose(" %d iterations %12.8lf gap\n", inner_count, e_upper -
e_lower);
e_lower = e_upper;
e_upper = 2 * GREEDY_BIG_E;
}
OUT:
log_debug(" %6d IPs (%.2lfms per call, %.2lfs total)\n", cplex_count,
cplex_time * 1000.0 / cplex_count, cplex_time);
for(int i = 0; i < nrays; i++)
abort_if(DOUBLE_iszero(bounds[i]), "bounds should be positive");
if_verbose_level
{
time_printf("Bounds:\n");
for (int k = 0; k < nrays; k++)
time_printf(" %12.8lf %12.8lf\n", k, bounds[k], 1 / bounds[k]);
}
CLEANUP:
if (x) free(x);
return rval;
}

1331
infinity/library/src/greedy-nd.c
View File

File diff suppressed because it is too large Load Diff

478
infinity/library/src/greedy.c
Unescape View File

@ -0,0 +1,478 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <stdlib.h>
#include <multirow/cg.h>
#include <multirow/double.h>
#include <multirow/util.h>
#include <infinity/greedy.h>
#include <infinity/greedy-bsearch.h>
#include <infinity/greedy-2d.h>
#include <infinity/greedy-nd.h>
struct SortPair
{
int index;
void *data;
};
static int _qsort_cmp_rays_angle(const void *p1,
const void *p2)
{
double *r1 = (double *) (((struct SortPair *) p1)->data);
double *r2 = (double *) (((struct SortPair *) p2)->data);
return sign(atan2(r1[0], r1[1]) - atan2(r2[0], r2[1]));
}
static int sort_rays_angle(double *rays,
int nrays,
double *beta)
{
int rval = 0;
int *map = 0;
double *rays_copy = 0;
double *beta_copy = 0;
struct SortPair *pairs = 0;
pairs = (struct SortPair *) malloc(nrays * sizeof(struct SortPair));
rays_copy = (double *) malloc(2 * nrays * sizeof(double));
beta_copy = (double*) malloc(nrays * sizeof(double));
abort_if(!pairs, "could not allocate pairs");
abort_if(!rays_copy, "could not allocate rays_copy");
abort_if(!beta_copy, "could not allocate beta_copy");
memcpy(rays_copy, rays, 2 * nrays * sizeof(double));
memcpy(beta_copy, beta, nrays * sizeof(double));
for (int i = 0; i < nrays; i++)
{
pairs[i].index = i;
pairs[i].data = &rays[2 * i];
}
qsort(pairs, nrays, sizeof(struct SortPair), _qsort_cmp_rays_angle);
for (int i = 0; i < nrays; i++)
{
beta[i] = beta_copy[pairs[i].index];
memcpy(&rays[2 * i], &rays_copy[2 * pairs[i].index], 2 *
sizeof(double));
}
CLEANUP:
if (pairs) free(pairs);
if (rays_copy) free(rays_copy);
if (beta_copy) free(beta_copy);
return rval;
}
static int scale_rays(double *rays,
int nrays,
double *scale)
{
int rval = 0;
for (int i = 0; i < nrays; i++)
{
rays[2*i] *= scale[i];
rays[2*i+1] *= scale[i];
}
CLEANUP:
return rval;
}
static void print_row(const struct Row *row)
{
time_printf("Row:\n");
for (int i = 0; i < row->nz; i++)
time_printf(" %.4lfx%d\n", row->pi[i], row->indices[i]);
time_printf(" <= %.4lf [%d]\n", row->pi_zero, row->head);
}
static void print_rays(const int *map,
const int *indices,
const double *f,
const double *rays,
const double *scale,
int nrays,
int nz,
int nrows)
{
time_printf("Mapping:\n");
for (int i = 0; i < nz; i++)
time_printf(" %4d: %4d x%-4d (scale=%.4lf)\n", i, map[i], indices[i],
scale[i]);
time_printf("Origin:\n");
for (int i = 0; i < nrows; i++)
time_printf(" %20.12lf\n", f[i]);
time_printf("Rays:\n");
for (int i = 0; i < nrays; i++)
{
time_printf(" ");
for (int j = 0; j < nrows; j++)
printf("%20.12lf ", rays[i * nrows + j]);
printf("angle=%.4lf ", atan2(rays[i * nrows], rays[i * nrows + 1]));
printf("norm=%.4lf ", fabs(rays[i * nrows]) + fabs(rays[i * nrows + 1]));
printf("[ ");
for (int j = 0; j < nz; j++)
if (map[j] == i) printf("%d ", indices[j]);
printf("]\n");
}
}
static void print_cut(const struct Row *cut)
{
time_printf("Generated cut:\n");
for (int i = 0; i < cut->nz; i++)
time_printf(" %.4lfx%d\n", cut->pi[i], cut->indices[i]);
time_printf(" <= %.4lf\n", cut->pi_zero);
}
int GREEDY_write_sage_file(int nrows,
int nrays,
const double *f,
const double *rays,
const double *beta,
const char *filename)
{
int rval = 0;
FILE *fsage = fopen(filename, "w");
abort_iff(!fsage, "could not open %s", filename);
fprintf(fsage, "f=vector([");
for (int i = 0; i < nrows; i++)
fprintf(fsage, "%.20lf,", f[i]);
fprintf(fsage, "])\n");
fprintf(fsage, "R=matrix([\n");
for (int i = 0; i < nrays; i++)
{
fprintf(fsage, " [");
for (int j = 0; j < nrows; j++)
{
fprintf(fsage, "%.20lf,", rays[i * nrows + j]);
}
fprintf(fsage, "],\n");
}
fprintf(fsage, "])\n");
fprintf(fsage, "pi=vector([\n");
for (int k = 0; k < nrays; k++)
fprintf(fsage, " %.12lf,\n", 1 / beta[k]);
fprintf(fsage, "])\n");
CLEANUP:
fclose(fsage);
return rval;
}
static int create_cut(const int nrows,
const char *column_types,
const int nrays,
const double *f,
const int lfree_nrays,
const double *lfree_rays,
const double *beta,
const double *extracted_rays,
const int nvars,
const int *variable_to_ray,
const int *indices,
const double *scale,
struct Row *cut)
{
int rval = 0;
cut->nz = nvars;
cut->pi = (double *) malloc(nvars * sizeof(double));
cut->indices = (int *) malloc(nvars * sizeof(int));
abort_if(!cut->pi, "could not allocate cut->pi");
abort_if(!cut->indices, "could not allocate cut->indices");
struct LP lp;
rval = LP_open(&lp);
abort_if(rval, "LP_open failed");
rval = GREEDY_create_psi_lp(nrows, lfree_nrays, f, lfree_rays, beta, &lp);
abort_if(rval, "create_psi_lp failed");
for (int i = 0; i < nvars; i++)
{
double value;
const double *q = &extracted_rays[variable_to_ray[i] * nrows];
if(ENABLE_LIFTING && column_types[indices[i]] == MILP_INTEGER)
{
rval = GREEDY_ND_pi(nrows, lfree_nrays, f, lfree_rays, beta, q,
scale[i], &lp, &value);
abort_if(rval, "GREEDY_ND_pi failed");
}
else
{
rval = GREEDY_ND_psi(nrows, lfree_nrays, f, lfree_rays, beta, q,
scale[i], &lp, &value);
abort_if(rval, "GREEDY_ND_psi failed");
}
log_verbose(" psi[%4d] = %20.12lf %d\n", indices[i], value);
value *= 1.0001;
value = DOUBLE_max(value, 0.0001);
cut->indices[i] = indices[i];
cut->pi[i] = - value;
}
cut->pi_zero = -1.0;
CLEANUP:
LP_free(&lp);
return rval;
}
#ifndef TEST_SOURCE
int GREEDY_generate_cut(int nrows,
struct Row **rows,
const char *column_types,
struct Row *cut)
{
int rval = 0;
double *f = 0;
long max_nrays = 0;
f = (double *) malloc(nrows * sizeof(double));
abort_if(!f, "could not allocate f");
for (int i = 0; i < nrows; i++)
{
f[i] = frac(rows[i]->pi_zero);
if (DOUBLE_eq(f[i], 1.0)) f[i] = 0.0;
max_nrays += rows[i]->nz;
}
int nz;
int nrays;
int *variable_to_ray = 0;
int *indices = 0;
double *extracted_rays = 0;
double *scale = 0;
double *beta = 0;
double *scaled_rays = 0;
int lfree_nrays = 0;
double *lfree_rays = 0;
variable_to_ray = (int *) malloc(max_nrays * sizeof(int));
indices = (int *) malloc(max_nrays * sizeof(int));
scale = (double *) malloc(max_nrays * sizeof(double));
extracted_rays = (double *) malloc(nrows * max_nrays * sizeof(double));
abort_if(!variable_to_ray, "could not allocate variable_to_ray");
abort_if(!indices, "could not allocate indices");
abort_if(!scale, "could not allocate scale");
abort_if(!extracted_rays, "could not allocate extracted_rays");
lfree_rays = (double*) malloc(nrows * (max_nrays + 100) * sizeof(double));
abort_if(!lfree_rays, "could not allocate lfree_rays");
rval = CG_extract_rays_from_rows(nrows, rows, extracted_rays, &nrays,
variable_to_ray, scale, indices, &nz);
abort_if(rval, "CG_extract_rays_from_rows failed");
for (double norm_cutoff = 0.00; norm_cutoff <= 5.0; norm_cutoff+= 0.1)
{
lfree_nrays = 0;
for (int i = 0; i < nrays; i++)
{
int keep = 1;
double *r = &extracted_rays[nrows * i];
for (int j = 0; j < lfree_nrays; j++)
{
double *q = &extracted_rays[nrows * j];
double norm = 0;
for(int k = 0; k < nrows; k++)
norm += fabs(r[k] - q[k]);
if(norm <= norm_cutoff)
{
keep = 0;
break;
}
}
if(keep)
{
memcpy(&lfree_rays[nrows * lfree_nrays], r, nrows * sizeof(double));
lfree_nrays++;
}
}
log_debug(" norm_cutoff=%8.2lf nrays=%8d\n", norm_cutoff, lfree_nrays);
if(lfree_nrays < MAX_N_RAYS) break;
}
if(ENABLE_LIFTING)
{
abort_if(nrows > 3, "not implemented");
int n_extra_rays;
double extra_rays[100];
if(nrows == 2)
{
extra_rays[0] = 0.0001;
extra_rays[1] = 0.0000;
extra_rays[2] = -0.0001;
extra_rays[3] = 0.0001;
extra_rays[4] = -0.0001;
extra_rays[5] = -0.0001;
n_extra_rays = 3;
}
else if(nrows == 3)
{
extra_rays[0] = 0.0000;
extra_rays[1] = 0.0000;
extra_rays[2] = 0.0001;
extra_rays[3] = 0.0001;
extra_rays[4] = 0.0000;
extra_rays[5] = -0.0001;
extra_rays[6] = -0.0001;
extra_rays[7] = 0.0000;
extra_rays[8] = -0.0001;
extra_rays[ 9] = -0.0001;
extra_rays[10] = -0.0001;
extra_rays[11] = -0.0001;
n_extra_rays = 4;
}
for(int i = 0; i < n_extra_rays; i++)
{
double *r = &extra_rays[nrows * i];
double scale;
int found, index;
rval = CG_find_ray(nrows, lfree_rays, lfree_nrays, r, &found,
&scale, &index);
abort_if(rval, "CG_find_ray failed");
if(!found) {
memcpy(&lfree_rays[nrows * lfree_nrays], r, nrows *
sizeof(double));
lfree_nrays++;
}
}
}
if(lfree_nrays < 3)
{
rval = ERR_NO_CUT;
cut->pi = 0;
cut->indices = 0;
goto CLEANUP;
}
log_debug("Extracted %d rays\n", lfree_nrays);
if_verbose_level
{
print_rays(variable_to_ray, indices, f, extracted_rays, scale, nrays,
nz, nrows);
}
beta = (double *) malloc((lfree_nrays) * sizeof(double));
abort_if(!beta, "could not allocate beta");
log_verbose("Computing lattice-free set...\n");
if(nrows == 2)
{
rval = sort_rays_angle(lfree_rays, lfree_nrays, beta);
abort_if(rval, "sort_rays_angle failed");
rval = GREEDY_2D_generate_cut(lfree_rays, lfree_nrays, f, beta);
if(rval)
{
rval = ERR_NO_CUT;
goto CLEANUP;
}
}
else
{
rval = GREEDY_ND_generate_cut(nrows, lfree_nrays, f, lfree_rays, beta);
if(rval)
{
rval = ERR_NO_CUT;
goto CLEANUP;
}
}
if(DUMP_CUT)
{
char filename[100];
sprintf(filename, "cut-%03d.sage", DUMP_CUT_N++);
time_printf("Writing %s...\n", filename);
rval = GREEDY_write_sage_file(nrows, lfree_nrays, f, lfree_rays, beta,
filename);
abort_if(rval, "GREEDY_write_sage_file failed");
}
rval = create_cut(nrows, column_types, nrays, f, lfree_nrays, lfree_rays,
beta, extracted_rays, nz, variable_to_ray, indices, scale, cut);
abort_if(rval, "create_cut failed");
if_verbose_level print_cut(cut);
CLEANUP:
if (f) free(f);
if (variable_to_ray) free(variable_to_ray);
if (beta) free(beta);
if (indices) free(indices);
if (scale) free(scale);
if (scaled_rays) free(scaled_rays);
if (lfree_rays) free(lfree_rays);
return rval;
}
#endif // TEST_SOURCE

616
infinity/library/tests/greedy-2d-test.cpp
Unescape View File

@ -0,0 +1,616 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <gtest/gtest.h>
#define TEST_SOURCE
extern "C" {
#include <multirow/util.h>
#include <infinity/greedy-2d.h>
#include "../src/greedy-2d.c"
}
#define BOUNDS_EPSILON 0.01
TEST(Greedy2DTest, test_generate_cut_1)
{
int rval = 0;
double bounds[100];
double f[] = {1 / 4.0, 3 / 4.0};
double rays[] = {-2 / 5.0, 5 / 7.0, 0.0, 1.0, 1.0, 1.0, 4 / 5.0, -2 / 3.0,
-1.0, 0.0};
rval = GREEDY_2D_generate_cut(rays, 5, f, bounds);
abort_if(rval, "GREEDY_2D_generate_cut failed");
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);
CLEANUP:
if (rval) FAIL();
}
TEST(Greedy2DTest, test_generate_cut_2)
{
int rval = 0;
double bounds[100];
double f[] = {1 / 2.0, 1 / 2.0};
double rays[] = {
-1.0, -1.0,
-1.0, 1.0,
1.0, 1.0,
1.0, 0.0,
1.0, -1.0
};
rval = GREEDY_2D_generate_cut(rays, 5, f, bounds);
abort_if(rval, "GREEDY_2D_generate_cut 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);
CLEANUP:
if (rval) FAIL();
}
TEST(Greedy2DTest, test_generate_cut_3)
{
int rval = 0;
double bounds[100];
double f[] = {5 / 22.0, 0.0};
double rays[] = {-1 / 22.0, 0.0, 0.0, 1 / 18.0, 1 / 22.0, 0.0};
rval = GREEDY_2D_generate_cut(rays, 3, f, bounds);
abort_if(rval, "GREEDY_2D_generate_cut failed");
EXPECT_NEAR(5.0, bounds[0], BOUNDS_EPSILON);
EXPECT_NEAR(17.0, bounds[2], BOUNDS_EPSILON);
EXPECT_EQ(GREEDY_BIG_E, bounds[1]);
CLEANUP:
if (rval) FAIL();
}
TEST(Greedy2DTest, scale_to_chull_test)
{
int rval = 0;
double rays[] = {
0, 1,
1, 1,
0.5, 0.25,
1, 0,
-1, -1,
-0.25, 0
};
double scale[6];
int nrays = 6;
rval = scale_to_chull(rays, nrays, scale);
abort_if(rval, "scale_to_chull failed");
EXPECT_NEAR(rays[0], 0.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[1], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[2], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[3], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[4], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[5], 0.5, BOUNDS_EPSILON);
EXPECT_NEAR(rays[6], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[7], 0.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[8],-1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[9],-1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[10],-0.5, BOUNDS_EPSILON);
EXPECT_NEAR(rays[11], 0.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[0], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[1], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[2], 2.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[3], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[4], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[5], 2.0, BOUNDS_EPSILON);
CLEANUP:
if(rval) FAIL();
}
TEST(Greedy2DTest, scale_to_chull_test_2)
{
int rval = 0;
double rays[] = {
1, 1,
0.5, 0.25,
1, 0,
};
double scale[3];
int nrays = 3;
rval = scale_to_chull(rays, nrays, scale);
abort_if(rval, "scale_to_chull failed");
EXPECT_NEAR(rays[0], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[1], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[2], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[3], 0.5, BOUNDS_EPSILON);
EXPECT_NEAR(rays[4], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(rays[5], 0.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[0], 1.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[1], 2.0, BOUNDS_EPSILON);
EXPECT_NEAR(scale[2], 1.0, BOUNDS_EPSILON);
CLEANUP:
if(rval) FAIL();
}
TEST(Greedy2DTest, find_containing_cone_test)
{
int rval = 0;
double rays[] = {-1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, -1.0, -1.0, -1.0};
double p1[] = {-1.0, 1.0};
double p2[] = {0.0, 1.0};
double p3[] = {1.0, 0.0};
double p4[] = {-1.0, 0.0};
int index1, index2;
double lambda1, lambda2;
rval = find_containing_cone(rays, 5, p1, &index1, &index2, &lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
EXPECT_EQ(0, index1);
EXPECT_EQ(1, index2);
rval = find_containing_cone(rays, 5, p2, &index1, &index2, &lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
EXPECT_EQ(0, index1);
EXPECT_EQ(1, index2);
rval = find_containing_cone(rays, 5, p3, &index1, &index2, &lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
EXPECT_EQ(1, index1);
EXPECT_EQ(2, index2);
rval = find_containing_cone(rays, 5, p4, &index1, &index2, &lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
EXPECT_EQ(4, index1);
EXPECT_EQ(0, index2);
CLEANUP:
if(rval) FAIL();
}
TEST(Greedy2DTest, find_containing_cone_test_2)
{
int rval = 0;
double rays[] = {0.0, 0.5, 1.0, 0.0, -0.5, 0.5};
double p[] = {0.5, -0.5};
int index1, index2;
double lambda1, lambda2;
rval = find_containing_cone(rays, 3, p, &index1, &index2, &lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
EXPECT_EQ(-1, index1);
EXPECT_EQ(-1, index2);
CLEANUP:
if(rval) FAIL();
}
TEST(Greedy2DTest, find_containing_cone_test_3)
{
int rval = 0;
double rays[] = {
-1, -0.1,
1, -0.1,
1, -0.5,
-1, -0.5,
};
double p[] = { 0, -2 };
int index1, index2;
double lambda1, lambda2;
rval = find_containing_cone(rays, 4, p, &index1, &index2, &lambda1, &lambda2);
abort_if(rval, "find_containing_cone failed");
EXPECT_EQ(2, index1);
EXPECT_EQ(3, index2);
CLEANUP:
if(rval) FAIL();
}
//TEST(Greedy2DTest, test_generate_cut_4)
//{
// int rval = 0;
// double bounds[100];
// double f[] = {5 / 22.0, 10 / 19.0};
// double rays[] = {0, -1 / 38.0, -1 / 22.0, -1 / 38.0, 0, 1 / 38.0, -1 / 22.0,
// 0, 1 / 22.0, 0, 1 / 22.0, 1 / 38.0};
//
// rval = GREEDY_2D_generate_cut(rays, 6, f, bounds);
// abort_if(rval, "GREEDY_2D_generate_cut failed");
//
// EXPECT_NEAR(20.0, bounds[0], BOUNDS_EPSILON);
// EXPECT_NEAR(20.0, bounds[1], BOUNDS_EPSILON);
// EXPECT_NEAR(18.0, bounds[2], BOUNDS_EPSILON);
// EXPECT_NEAR(18.0, bounds[5], BOUNDS_EPSILON);
// EXPECT_EQ(GREEDY_BIG_E, bounds[3]);
// EXPECT_EQ(GREEDY_BIG_E, bounds[4]);
//
// CLEANUP:
// if (rval) FAIL();
//}
//
//TEST(Greedy2DTest, test_generate_cut_5)
//{
// int rval = 0;
// double bounds[100];
// double f[] = {0.22727272727272729291, 0.52631578947368418131};
// double rays[] = {0.00000000000000000000, -0.02631578947368420907,
// -0.04545454545454545581, -0.02631578947368420907,
// 0.00000000000000000000, 0.02631578947368420907,
// -0.04545454545454545581, 0.00000000000000000000,
// 0.04545454545454545581, 0.00000000000000000000,
// 0.04545454545454545581, 0.02631578947368420907};
//
// rval = GREEDY_2D_generate_cut(rays, 6, f, bounds);
// abort_if(rval, "GREEDY_2D_generate_cut failed");
//
// EXPECT_NEAR(20.0, bounds[0], BOUNDS_EPSILON);
// EXPECT_NEAR(20.0, bounds[1], BOUNDS_EPSILON);
// EXPECT_NEAR(18.0, bounds[2], BOUNDS_EPSILON);
// EXPECT_NEAR(18.0, bounds[5], BOUNDS_EPSILON);
// EXPECT_EQ(GREEDY_BIG_E, bounds[3]);
// EXPECT_EQ(GREEDY_BIG_E, bounds[4]);
//
// CLEANUP:
// if (rval) FAIL();
//}
//
//
//
//TEST(Greedy2DTest, get_peak_ray_test_1)
//{
// int rval = 0;
//
// double rays[] = {-2.0, 0.0, -1.0, 1.0, 1.0, 1.0, 2.0, 0.0, 1.0, -1.0, -1.0, -1.0 };
// int nrays = 6;
//
// double normal1[] = { 0.0, 1.0 };
// double normal2[] = { 1.0, 1.0 };
// double normal3[] = { 0.0, -1.0 };
// double normal4[] = {-1.0, -1.0 };
// int index;
//
// rval = get_peak_ray(rays, nrays, normal1, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_TRUE(index == 1 || index == 2);
//
// rval = get_peak_ray(rays, nrays, normal3, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_TRUE(index == 4 || index == 5);
//
// rval = get_peak_ray(rays, nrays, normal2, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(2, index);
//
// rval = get_peak_ray(rays, nrays, normal4, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(5, index);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(Greedy2DTest, get_peak_ray_test_2)
//{
// int rval = 0;
//
// double rays[] = { -1e100, 0, 0, 1, 1e100, 0, 1, -1, -1,-1 };
// int nrays = 5;
//
// double normal1[] = { 0.0, 1.0 };
// double normal2[] = { 1.0, 1.0 };
// double normal3[] = { 1.0, 0.0 };
// double normal4[] = { 0.0, -1.0 };
// double normal5[] = { -1.0, -1.0 };
// double normal6[] = { -1.0, 0.0 };
// int index;
//
// rval = get_peak_ray(rays, nrays, normal1, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(1, index);
//
// rval = get_peak_ray(rays, nrays, normal4, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(3, index);
//
// rval = get_peak_ray(rays, nrays, normal2, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(2, index);
//
// rval = get_peak_ray(rays, nrays, normal5, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(0, index);
//
// rval = get_peak_ray(rays, nrays, normal3, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(2, index);
//
// rval = get_peak_ray(rays, nrays, normal6, &index);
// abort_if(rval, "get_peak_ray failed");
// EXPECT_EQ(0, index);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(Greedy2DTest, get_peak_ray_test_3)
//{
// int rval = 0;
//
// double rays[] = { -1, 1, 0, 1, 1, 1};
// int nrays = 3;
//
// double normal1[] = { 0.0, -1.0 };
// int index;
//
// rval = get_peak_ray(rays, nrays, normal1, &index);
// EXPECT_EQ(ERR_NOT_FOUND, rval);
// rval = 0;
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(Greedy2DTest, nearest_lattice_point_test_1)
//{
// int rval = 0;
//
// double z0a[] = {100, 1};
// double z0b[] = {-100, 1};
// double g[] = {-1, 0};
//
// double x1[] = {-1.34, 1};
// double x2[] = {-3, 1};
//
// double z1[2];
// double z2[2];
//
// rval = nearest_lattice_points(g, z0a, x1, z1, z2);
// abort_if(rval, "nearest_lattice_points failed");
// EXPECT_DOUBLE_EQ(-1, z1[0]);
// EXPECT_DOUBLE_EQ(1, z1[1]);
// EXPECT_DOUBLE_EQ(-2, z2[0]);
// EXPECT_DOUBLE_EQ(1, z2[1]);
//
// rval = nearest_lattice_points(g, z0b, x1, z1, z2);
// abort_if(rval, "nearest_lattice_points failed");
// EXPECT_DOUBLE_EQ(-1, z1[0]);
// EXPECT_DOUBLE_EQ(1, z1[1]);
// EXPECT_DOUBLE_EQ(-2, z2[0]);
// EXPECT_DOUBLE_EQ(1, z2[1]);
//
// rval = nearest_lattice_points(g, z0a, x2, z1, z2);
// abort_if(rval, "nearest_lattice_points failed");
// EXPECT_DOUBLE_EQ(-3, z1[0]);
// EXPECT_DOUBLE_EQ(1, z1[1]);
// EXPECT_DOUBLE_EQ(-3, z2[0]);
// EXPECT_DOUBLE_EQ(1, z2[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(Greedy2DTest, nearest_lattice_point_test_2)
//{
// int rval = 0;
//
// double z0[] = {1, 1};
// double g[] = {1, 2};
// double x[] = {3.35, 5.7};
//
// double z1[2];
// double z2[2];
//
// rval = nearest_lattice_points(g, z0, x, z1, z2);
// abort_if(rval, "nearest_lattice_points failed");
// EXPECT_DOUBLE_EQ(3, z1[0]);
// EXPECT_DOUBLE_EQ(5, z1[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(Greedy2DTest, find_normal_test)
//{
// int rval = 0;
//
// double d1[] = {1.0, 0.0};
// double d2[] = {1.0, 1.0};
//
// double x1[] = {1.0, 1.0};
// double x2[] = {1.0, -1.0};
//
// double n[2];
//
// rval = find_normal(d1, x1, n);
// abort_if(rval, "find_normal failed");
// EXPECT_DOUBLE_EQ(0, n[0]);
// EXPECT_DOUBLE_EQ(1, n[1]);
//
// rval = find_normal(d1, x2, n);
// abort_if(rval, "find_normal failed");
// EXPECT_DOUBLE_EQ(0, n[0]);
// EXPECT_DOUBLE_EQ(-1, n[1]);
//
// rval = find_normal(d2, x1, n);
// abort_if(rval, "find_normal failed");
// EXPECT_DOUBLE_EQ(0, n[0]);
// EXPECT_DOUBLE_EQ(0, n[1]);
//
// rval = find_normal(d2, x2, n);
// abort_if(rval, "find_normal failed");
// EXPECT_DOUBLE_EQ(1, n[0]);
// EXPECT_DOUBLE_EQ(-1, n[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(Greedy2DTest, check_rays_parallel)
//{
// double r1[] = {1.0, 2.0};
// double r2[] = {2.0, 4.0};
// double r3[] = {-2.0, -4.0};
// double r4[] = {0.0, 0.0};
//
// int match;
// double scale;
//
// check_rays_parallel(r1, r2, &match, &scale);
// EXPECT_TRUE(match);
// EXPECT_DOUBLE_EQ(scale, 0.5);
//
// check_rays_parallel(r1, r1, &match, &scale);
// EXPECT_TRUE(match);
// EXPECT_DOUBLE_EQ(scale, 1.0);
//
// check_rays_parallel(r2, r3, &match, &scale);
// EXPECT_FALSE(match);
//
// check_rays_parallel(r1, r4, &match, &scale);
// EXPECT_FALSE(match);
//}
//
//TEST(Greedy2DTest, find_ray)
//{
// double rays[] = {1.0, 2.0, -1.0, 0.0, -5.0, -5.0};
// int nrays = 3;
//
// double r1[] = {1.0, 2.0};
// double r2[] = {-2.0, 0.0};
// double r3[] = {-1.0, -1.0};
// double r4[] = {7.0, 1.0};
// double r5[] = {-1.0, -2.0};
//
// int found, index;
// double scale;
//
// find_ray(rays, nrays, r1, &found, &scale, &index);
// EXPECT_TRUE(found);
// EXPECT_EQ(index, 0);
// EXPECT_DOUBLE_EQ(scale, 1.0);
//
// find_ray(rays, nrays, r2, &found, &scale, &index);
// EXPECT_TRUE(found);
// EXPECT_EQ(index, 1);
// EXPECT_DOUBLE_EQ(scale, 2.0);
//
// find_ray(rays, nrays, r3, &found, &scale, &index);
// EXPECT_TRUE(found);
// EXPECT_EQ(index, 2);
// EXPECT_DOUBLE_EQ(scale, 1 / 5.0);
//
// find_ray(rays, nrays, r4, &found, &scale, &index);
// EXPECT_FALSE(found);
//
// find_ray(rays, nrays, r5, &found, &scale, &index);
// EXPECT_FALSE(found);
//}
//
//TEST(Greedy2DTest, extract_rays_from_two_sparse_rows_test)
//{
// int rval = 0;
//
// char ctypes[] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
// double pi1[] = {-3.0, 1.0, -2.0, -1.5, -3.0};
// int indices1[] = {2, 3, 6, 7, 9};
//
// double pi2[] = {1.0, -2.0, -3.0, -1.0};
// int indices2[] = {0, 2, 5, 7};
//
// struct Row row1 = {.nz = 5, .head = 3, .pi_zero = 0, .pi = pi1, .indices = indices1};
// struct Row row2 = {.nz = 4, .head = 0, .pi_zero = 0, .pi = pi2, .indices = indices2};
//
// int nrays;
// int nz;
//
// double rays[100];
// int indices[100];
// double ray_scale[100];
// int variable_to_ray[100];
//
// rval = extract_rays_from_two_sparse_rows(&row1, &row2, ctypes, rays, &nrays,
// variable_to_ray, ray_scale, indices, &nz);
// abort_if(rval, "extract_rays_from_two_sparse_rows failed");
//
// EXPECT_EQ(nrays, 3);
// EXPECT_EQ(nz, 5);
//
// EXPECT_DOUBLE_EQ(rays[0], 3.0);
// EXPECT_DOUBLE_EQ(rays[1], 2.0);
// EXPECT_DOUBLE_EQ(rays[2], 0.0);
// EXPECT_DOUBLE_EQ(rays[3], 3.0);
// EXPECT_DOUBLE_EQ(rays[4], 2.0);
// EXPECT_DOUBLE_EQ(rays[5], 0.0);
//
// EXPECT_EQ(indices[0], 2);
// EXPECT_EQ(indices[1], 5);
// EXPECT_EQ(indices[2], 6);
// EXPECT_EQ(indices[3], 7);
// EXPECT_EQ(indices[4], 9);
//
// EXPECT_EQ(variable_to_ray[0], 0);
// EXPECT_EQ(variable_to_ray[1], 1);
// EXPECT_EQ(variable_to_ray[2], 2);
// EXPECT_EQ(variable_to_ray[3], 0);
// EXPECT_EQ(variable_to_ray[4], 2);
//
// EXPECT_DOUBLE_EQ(ray_scale[0], 1.0);
// EXPECT_DOUBLE_EQ(ray_scale[1], 1.0);
// EXPECT_DOUBLE_EQ(ray_scale[2], 1.0);
// EXPECT_DOUBLE_EQ(ray_scale[3], 0.5);
// EXPECT_DOUBLE_EQ(ray_scale[4], 1.5);
//
// CLEANUP:
// if (rval) FAIL();
//}
//

457
infinity/library/tests/greedy-nd-test.cpp
Unescape View File

@ -0,0 +1,457 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <gtest/gtest.h>
extern "C" {
#include <math.h>
#include <multirow/lp.h>
#include <multirow/util.h>
#include <infinity/greedy-nd.h>
}
int ENABLE_LIFTING = 0;
int MIN_N_ROWS = 2;
int MAX_N_ROWS = 2;
int DUMP_CUT = 0;
int DUMP_CUT_N = 0;
TEST(GreedyNDTest, find_violated_cone_test)
{
int rval = 0;
int nrows = 2;
int nrays = 4;
double f[] = { 0.5, 0.5 };
double rays[] =
{
-1.0, 1.0,
1.0, 1.0,
1.0, -1.0,
-1.0, -1.0
};
double x[] = { 1.0, 0.5 };
double beta[] = { INFINITY, INFINITY, INFINITY, INFINITY };
double sbar[nrays];
int rx[nrays];
int found;
rval = GREEDY_ND_find_violated_cone(nrows, nrays, f, rays, x, beta,
1.0, rx, sbar, &found);
abort_if(rval, "GREEDY_ND_find_violated_cone failed");
EXPECT_TRUE(found);
EXPECT_FALSE(rx[0]);
EXPECT_TRUE(rx[1]);
EXPECT_TRUE(rx[2]);
EXPECT_FALSE(rx[3]);
rval = GREEDY_ND_find_violated_cone(nrows, nrays, f, rays, x, beta,
0.5, rx, sbar, &found);
abort_if(rval, "GREEDY_ND_find_violated_cone failed");
EXPECT_FALSE(found);
CLEANUP:
if(rval) FAIL();
}
TEST(GreedyNDTest, find_tight_rays_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 beta[] = { INFINITY, INFINITY, INFINITY, INFINITY, INFINITY,
INFINITY };
double epsilon = 0.5;
double x[] = { 1.0, 1.0 };
int tx[6];
rval = GREEDY_ND_find_tight_rays(2, 6, f, rays, x, beta, epsilon, tx);
abort_if(rval, "GREEDY_ND_find_tight_rays failed");
EXPECT_TRUE(tx[0]);
EXPECT_FALSE(tx[1]);
EXPECT_FALSE(tx[2]);
EXPECT_FALSE(tx[3]);
EXPECT_FALSE(tx[4]);
EXPECT_FALSE(tx[5]);
CLEANUP:
if(rval) FAIL();
}
TEST(GreedyNDTest, find_tight_rays_test_2)
{
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 beta[] = { 0.5, 0.5, 0.5, 0.5, INFINITY, INFINITY };
double epsilon = 1.0;
double x[] = { 1.0, 1.0 };
int tx[6];
rval = GREEDY_ND_find_tight_rays(2, 6, f, rays, x, beta, epsilon, tx);
abort_if(rval, "GREEDY_ND_find_tight_rays failed");
EXPECT_TRUE(tx[0]);
EXPECT_FALSE(tx[1]);
EXPECT_FALSE(tx[2]);
EXPECT_FALSE(tx[3]);
EXPECT_TRUE(tx[4]);
EXPECT_TRUE(tx[5]);
CLEANUP:
if(rval) FAIL();
}
TEST(GreedyNDTest, cone_bound_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 beta[] = { INFINITY, INFINITY, INFINITY, INFINITY, INFINITY,
INFINITY };
double x[] = { 1.0, 1.0 };
int rx1[] = { 1, 0, 0, 0, 0, 0 };
int rx2[] = { 0, 0, 0, 0, 1, 1 };
double epsilon;
rval = GREEDY_ND_cone_bound(2, 6, f, rays, rx1, x, beta, &epsilon);
abort_if(rval, "GREEDY_ND_cone_bound failed");
EXPECT_NEAR(0.5, epsilon, 1e-6);
rval = GREEDY_ND_cone_bound(2, 6, f, rays, rx2, x, beta, &epsilon);
abort_if(rval, "GREEDY_ND_cone_bound failed");
EXPECT_NEAR(1.0, epsilon, 1e-6);
CLEANUP:
if(rval) FAIL();
}
TEST(GreedyNDTest, cone_bound_test_2)
{
int rval = 0;
double f[] = { 0.0, 0.0 };
double rays[] =
{
0.0, 1.0,
1.0, 0.0
};
int rx[] = { 1, 1 };
double beta1[] = { 100, 100 };
double beta2[] = { 0.5, 100 };
double beta3[] = { 0.5, 1.0 };
double x1[] = { 0.5, 0.5 };
double x2[] = { 0.5, 0.25 };
double epsilon;
rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x1, beta1, &epsilon);
abort_if(rval, "GREEDY_ND_cone_bound failed");
EXPECT_NEAR(1.0, epsilon, 1e-6);
rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x1, beta2, &epsilon);
abort_if(rval, "GREEDY_ND_cone_bound failed");
EXPECT_EQ(INFINITY, epsilon);
rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x2, beta2, &epsilon);
abort_if(rval, "GREEDY_ND_cone_bound failed");
EXPECT_NEAR(1.0, epsilon, 1e-6);
rval = GREEDY_ND_cone_bound(2, 2, f, rays, rx, x2, beta3, &epsilon);
abort_if(rval, "GREEDY_ND_cone_bound failed");
EXPECT_EQ(INFINITY, epsilon);
CLEANUP:
if(rval) FAIL();
}
TEST(GreedyNDTest, bound_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 beta1[] = { INFINITY, INFINITY, INFINITY, INFINITY, INFINITY, INFINITY };
double beta2[] = { 0.5, 0.5, 0.5, 0.5, INFINITY, INFINITY };
double beta3[] = { 0.5, 0.5, 0.5, 0.5, 1.0, 1.0 };
double x[] = { 1.0, 1.0 };
double epsilon;
int tx[6];
rval = GREEDY_ND_bound(2, 6, f, rays, x, beta1, &epsilon, tx);
abort_if(rval, "GREEDY_ND_bound failed");
EXPECT_NEAR(epsilon, 0.5, 1e-6);
EXPECT_TRUE(tx[0]);
EXPECT_FALSE(tx[1]);
EXPECT_FALSE(tx[2]);
EXPECT_FALSE(tx[3]);
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");
EXPECT_NEAR(epsilon, 1.0, 1e-6);
EXPECT_TRUE(tx[0]);
EXPECT_FALSE(tx[1]);
EXPECT_FALSE(tx[2]);
EXPECT_FALSE(tx[3]);
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");
EXPECT_EQ(epsilon, INFINITY);
EXPECT_FALSE(tx[0]);
EXPECT_FALSE(tx[1]);
EXPECT_FALSE(tx[2]);
EXPECT_FALSE(tx[3]);
EXPECT_FALSE(tx[4]);
EXPECT_FALSE(tx[5]);
CLEANUP:
if(rval) FAIL();
}
TEST(GreedyNDTest, psi_test)
{
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 beta[] = { 0.5, 0.5, 0.5, 0.5, 1.0, 1.0 };
double q1[] = { 1.0, 1.0 };
double q2[] = { -2.0, 0.0 };
double value;
struct LP lp;
rval = LP_open(&lp);
abort_if(rval, "LP_open failed");
rval = GREEDY_create_psi_lp(2, 6, f, rays, beta, &lp);
abort_if(rval, "GREEDY_create_psi_lp failed");
rval = GREEDY_ND_psi(2, 6, f, rays, beta, q1, 1.0, &lp, &value);
abort_if(rval, "GREDDY_ND_psi failed");
EXPECT_NEAR(value, 2.0, 1e-6);
rval = GREEDY_ND_psi(2, 6, f, rays, beta, q2, 2.0, &lp, &value);
abort_if(rval, "GREDDY_ND_psi failed");
EXPECT_NEAR(value, 8.0, 1e-6);
CLEANUP:
LP_free(&lp);
if(rval) FAIL();
}
TEST(GreedyNDTest, psi_test_2)
{
int rval = 0;
double f[] = { 0.5, 0.5, 0.5 };
double rays[] =
{
-0.5, -0.5, -0.5,
-0.5, -0.5, 0.5,
-0.5, 0.5, -0.5,
-0.5, 0.5, 0.5,
0.5, -0.5, -0.5,
0.5, -0.5, 0.5,
0.5, 0.5, -0.5,
0.5, 0.5, 0.5,
};
double beta[] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
double q1[] = { 0.5, 0.5, 0.5 };
double q2[] = { 1, 0, 0 };
double value;
struct LP lp;
rval = LP_open(&lp);
abort_if(rval, "LP_open failed");
rval = GREEDY_create_psi_lp(3, 8, f, rays, beta, &lp);
abort_if(rval, "GREEDY_create_psi_lp failed");
rval = GREEDY_ND_psi(3, 8, f, rays, beta, q1, 1.0, &lp, &value);
abort_if(rval, "GREDDY_ND_psi failed");
EXPECT_NEAR(value, 1.0, 1e-6);
rval = GREEDY_ND_psi(3, 8, f, rays, beta, q2, 1.0, &lp, &value);
abort_if(rval, "GREDDY_ND_psi failed");
EXPECT_NEAR(value, 2.0, 1e-6);
CLEANUP:
if(rval) FAIL();
}
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 beta[6];
rval = GREEDY_ND_generate_cut(2, 6, f, rays, beta);
abort_if(rval, "GREEDY_ND_generate_cut failed");
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);
CLEANUP:
if(rval) FAIL();
}
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 beta[6];
rval = GREEDY_ND_generate_cut(3, 6, f, rays, beta);
abort_if(rval, "GREEDY_ND_generate_cut 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);
CLEANUP:
if(rval) FAIL();
}
TEST(GreedyNTest, scale_to_ahull_test)
{
int rval = 0;
double rays[] =
{
0.0, 0.0, 1.0,
0.0, 1.0, 0.0,
1.0, 0.0, 0.0,
-1.0, 0.0, 0.0
};
int rx[] = { 1, 1, 1 , 0 };
double beta[] = { INFINITY, INFINITY, INFINITY, INFINITY };
double epsilon = 1.0;
double d1[] = { 1.0, 1.0, 1.0 };
double d2[] = { 2.0, 2.0, 0.0 };
double d3[] = { -1.0, -1.0, -1.0 };
double alpha;
rval = GREEDY_ND_scale_to_ahull(3, 4, rays, rx, beta, epsilon, d1, &alpha);
abort_if(rval, "GREEDY_ND_scale_to_ahull failed");
EXPECT_DOUBLE_EQ(1 / 3.0, alpha);
rval = GREEDY_ND_scale_to_ahull(3, 4, rays, rx, beta, epsilon, d2, &alpha);
abort_if(rval, "GREEDY_ND_scale_to_ahull failed");
EXPECT_DOUBLE_EQ(0.25, alpha);
rval = GREEDY_ND_scale_to_ahull(3, 4, rays, rx, beta, epsilon, d3, &alpha);
abort_if(rval, "GREEDY_ND_scale_to_ahull failed");
EXPECT_EQ(INFINITY, alpha);
CLEANUP:
if(rval) FAIL();
}

2
lifting/benchmark/CMakeLists.txt
Unescape View File

@ -0,0 +1,2 @@
add_executable(lifting-benchmark.run src/main.c)
target_link_libraries (lifting-benchmark.run LINK_PUBLIC lifting_static multirow_static)

639
lifting/benchmark/src/main.c
Unescape View File

@ -0,0 +1,639 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define _XOPEN_SOURCE 500
#include <stdio.h>
#include <stdarg.h>
#include <getopt.h>
#include <unistd.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>
#include <stdlib.h>
#include <multirow/util.h>
#include <multirow/double.h>
#include <multirow/lfree2d.h>
#include <lifting/lifting.h>
#include <lifting/lifting-mip.h>
char LOG_FILENAME[1000] = {0};
char STATS_FILENAME[1000] = {0};
char SETS_FILENAME[1000] = {0};
char ANSWERS_FILENAME[1000] = {0};
unsigned int SEED = 0;
#define ALGORITHM_BOUND 0
#define ALGORITHM_NAIVE 1
#define ALGORITHM_MIP 2
int SELECT_NAIVE_ALGORITHM = 0;
int SELECT_BOUND_ALGORITHM = 0;
int SELECT_MIP_ALGORITHM = 0;
int ENABLE_PREPROCESSING = 0;
int ENABLE_SHEAR = 0;
int CHECK_ANSWERS = 0;
int WRITE_ANSWERS = 0;
int USE_FIXED_BOUNDS = 0;
int NAIVE_BIG_M = 0;
int N_SAMPLES_PER_SET = 10;
int N_ANSWERS = 0;
double ANSWERS[100000];
double PRE_M[100000];
double CENTER[100000];
FILE *LOG_FILE;
FILE *STATS_FILE;
FILE *ANSWERS_FILE;
int BOUNDING_BOX_PADDING = 5;
static const struct option options_tab[] =
{
{"help", no_argument, 0, 'h'},
{"sets", required_argument, 0, 'b'},
{"log", required_argument, 0, 'l'},
{"stats", required_argument, 0, 'o'},
{"seed", required_argument, 0, 's'},
{"fixed-bounds", required_argument, 0, 'f'},
{"naive", no_argument, 0, 'n'},
{"bound", no_argument, 0, 'u'},
{"preprocess", no_argument, 0, 'p'},
{"shear", no_argument, 0, 'e'},
{"write-answers", required_argument, 0, 'w'},
{"check-answers", required_argument, 0, 'c'},
{"samples", required_argument, 0, 'a'},
{"mip", no_argument, 0, 'm'},
{0, 0, 0, 0}
};
static void print_usage(char **argv)
{
printf("Usage: %s [OPTION]...\n", argv[0]);
printf("Performs trivial lifting using different algorithms.\n\n");
printf("Parameters:\n");
printf("%4s %-20s %s\n", "-l", "--log=FILE",
"write log to the specified file");
printf("%4s %-20s %s\n", "-b", "--sets=FILE",
"file containing lattice-free sets");
printf("%4s %-20s %s\n", "-o", "--stats=FILE",
"write statistics to the specified file");
printf("%4s %-20s %s\n", "-s", "--seed=NUM",
"random seed");
printf("%4s %-20s %s\n", "-f", "--fixed-bounds=NUM",
"user fixed bounds for naive algorithm");
printf("%4s %-20s %s\n", "-n", "--naive",
"select naive algorithm");
printf("%4s %-20s %s\n", "-m", "--mip",
"select MIP algorithm");
printf("%4s %-20s %s\n", "-u", "--bound",
"select bound algorithm");
printf("%4s %-20s %s\n", "-p", "--preprocess",
"enable pre-processing step in bound algorithm");
printf("%4s %-20s %s\n", "-e", "--shear",
"apply shear transformation to sets");
printf("%4s %-20s %s\n", "-w", "--write-answers=FILE",
"write computed coefficients to given file");
printf("%4s %-20s %s\n", "-c", "--check-answers=FILE",
"check computed coefficients against given file");
printf("%4s %-20s %s\n", "-a", "--samples=NUM",
"use specified number of samples per set");
}
void stats_printf(const char *fmt,
...)
{
if(STATS_FILE)
{
va_list args;
va_start(args, fmt);
vfprintf(STATS_FILE, fmt, args);
va_end(args);
fflush(STATS_FILE);
}
}
static int parse_args(int argc,
char **argv)
{
int rval = 0;
opterr = 0;
while (1)
{
int c = 0;
int option_index = 0;
c = getopt_long(argc, argv, "hb:k:s:f:o:nupew:c:a:m", options_tab,
&option_index);
if (c < 0) break;
switch (c)
{
case 'l':
strcpy(LOG_FILENAME, optarg);
break;
case 's':
{
int count = sscanf(optarg, "%u", &SEED);
abort_if(count != 1, "invalid seed");
break;
}
case 'a':
{
int count = sscanf(optarg, "%d", &N_SAMPLES_PER_SET);
abort_if(count != 1, "invalid number of samples");
abort_if(N_SAMPLES_PER_SET <= 0, "invalid number of samples");
break;
}
case 'f':
{
USE_FIXED_BOUNDS = 1;
int count = sscanf(optarg, "%d", &NAIVE_BIG_M);
abort_if(count != 1, "invalid fixed bound");
break;
}
case 'o':
strcpy(STATS_FILENAME, optarg);
break;
case 'b':
strcpy(SETS_FILENAME, optarg);
break;
case 'w':
strcpy(ANSWERS_FILENAME, optarg);
WRITE_ANSWERS = 1;
break;
case 'c':
strcpy(ANSWERS_FILENAME, optarg);
CHECK_ANSWERS = 1;
break;
case 'n':
SELECT_NAIVE_ALGORITHM = 1;
break;
case 'm':
SELECT_MIP_ALGORITHM = 1;
break;
case 'u':
SELECT_BOUND_ALGORITHM = 1;
break;
case 'p':
ENABLE_PREPROCESSING = 1;
break;
case 'e':
ENABLE_SHEAR = 1;
break;
case 'h':
print_usage(argv);
exit(0);
case ':':
fprintf(stderr, "%s: option '-%c' requires an argument\n",
argv[0], optopt);
rval = 1;
goto CLEANUP;
case '?':
default:
fprintf(stderr, "%s: option '-%c' is invalid\n", argv[0],
optopt);
rval = 1;
goto CLEANUP;
}
}
if ((strlen(SETS_FILENAME) == 0))
{
fprintf(stderr, "You must specify a file containing the lattice-free "
"sets.\n");
rval = 1;
}
if (SELECT_NAIVE_ALGORITHM + SELECT_BOUND_ALGORITHM + SELECT_MIP_ALGORITHM != 1)
{
fprintf(stderr, "You must select exactly one algorithm.\n");
rval = 1;
}
if (CHECK_ANSWERS + WRITE_ANSWERS > 1)
{
fprintf(stderr, "Cannot write and check answers at same time.\n");
rval = 1;
}
CLEANUP:
if (rval)
fprintf(stderr, "Try '%s --help' for more information.\n", argv[0]);
return rval;
}
void print_header(int argc,
char *const *argv)
{
char hostname[3000];
gethostname(hostname, 1024);
time_t now;
time(&now);
struct tm *ttm = localtime(&now);
time_printf("benchmark.run\n");
time_printf("%s %04d-%02d-%02d %02d:%02d\n", hostname, ttm->tm_year + 1900,
ttm->tm_mon + 1, ttm->tm_mday, ttm->tm_hour, ttm->tm_min);
time_printf("Compile-time parameters:\n");
time_printf(" EPSILON: %e\n", EPSILON);
char cmdline[5000] = {0};
for (int i = 0; i < argc; i++)
sprintf(cmdline + strlen(cmdline), "%s ", argv[i]);
time_printf("Command line arguments:\n");
time_printf(" %s\n", cmdline);
}
int benchmark_set_sample(int algorithm,
int set_idx,
const struct LFreeSet2D *set,
const double *rays,
const double *pre_m,
const double center,
int *lb,
int *ub,
int current_sample,
int *wrong_answer)
{
int rval = 0;
for (int i = 0; i < N_RAYS; i++)
{
double ray[2] = { rays[2 * i], rays[2 * i + 1] };
double value;
if(ENABLE_PREPROCESSING)
{
double r0 = ray[0], r1 = ray[1];
ray[0] = pre_m[0] * r0 + pre_m[2] * r1;
ray[1] = pre_m[1] * r0 + pre_m[3] * r1;
ray[0] = ray[0] - floor(set->f[0] + ray[0]);
ray[1] = ray[1] + floor(center + 0.5 - set->f[1] - ray[1]);
}
log_debug(" Ray %d (%.6lf,%.6lf)...\n", i, ray[0], ray[1]);
switch (algorithm)
{
case ALGORITHM_BOUND:
rval = LIFTING_2D_bound(set->n_halfspaces, set->halfspaces, ray,
&value);
abort_if(rval, "LIFTING_2D_bound failed");
break;
case ALGORITHM_NAIVE:
rval = LIFTING_2D_naive(set->n_halfspaces, set->halfspaces, ray,
lb, ub, &value);
abort_if(rval, "LIFTING_2D_naive failed");
break;
case ALGORITHM_MIP:
rval = LIFTING_2D_mip(set->n_halfspaces, set->halfspaces, ray,
&value);
abort_if(rval, "LIFTING_2D_mip failed");
break;
default:
abort_if(1, "Invalid algorithm");
}
if(current_sample == 0)
{
if(WRITE_ANSWERS)
{
abort_iff(!DOUBLE_geq(value, 0),
"value should be non-negative (%.8lf)", value);
fprintf(ANSWERS_FILE, "%d %d %.20lf\n", set_idx, i, value);
}
if(CHECK_ANSWERS)
{
int count;
int expected_set_idx, expected_i;
double expected_value;
while(1)
{
count = fscanf(ANSWERS_FILE, "%d %d %lf ",
&expected_set_idx, &expected_i, &expected_value);
abort_if(count != 3, "error reading answer");
if(set_idx == expected_set_idx && i == expected_i) break;
}
double delta = fabs(value - expected_value);
if(delta > 1e-3)
{
log_warn(" wrong answer (set=%d ray=%d answer=%.8lf"
" expected=%.8lf delta=%.8lf)\n", set_idx,
i, value, expected_value, delta);
*wrong_answer = 1;
LFREE_2D_print_set(set);
}
}
log_verbose(" %4d %12.8lf\n", j, value);
}
}
CLEANUP:
return rval;
}
int benchmark_set(int algorithm,
int set_idx,
const struct LFreeSet2D *set,
const double *rays,
const double *pre_m,
const double center,
int *wrong_answer)
{
int rval = 0;
int lb[2], ub[2];
if(algorithm == ALGORITHM_NAIVE)
{
if (USE_FIXED_BOUNDS)
{
ub[0] = ub[1] = NAIVE_BIG_M;
lb[0] = lb[1] = -NAIVE_BIG_M;
}
else
{
rval = LFREE_2D_get_bounding_box(set, lb, ub);
abort_if(rval, "LFREE_2D_get_bounding_box failed");
ub[0] += BOUNDING_BOX_PADDING;
ub[1] += BOUNDING_BOX_PADDING;
lb[0] -= BOUNDING_BOX_PADDING;
lb[0] -= BOUNDING_BOX_PADDING;
}
ub[0] = fmin(ub[0], MAX_BOX_SIZE);
ub[1] = fmin(ub[1], MAX_BOX_SIZE);
lb[0] = fmax(lb[0], -MAX_BOX_SIZE);
lb[1] = fmax(lb[1], -MAX_BOX_SIZE);
if(ub[0] == MAX_BOX_SIZE || ub[1] == MAX_BOX_SIZE || lb[0] ==
-MAX_BOX_SIZE || lb[1] == -MAX_BOX_SIZE)
{
log_info(" bounding box has been clipped\n");
}
}
for (int k = 0; k < N_SAMPLES_PER_SET; k ++)
{
rval = benchmark_set_sample(algorithm, set_idx, set, rays, pre_m,
center, lb, ub, k, wrong_answer);
abort_if(rval, "benchmark_set_sample failed");
}
CLEANUP:
return rval;
}
int benchmark(int n_sets, struct LFreeSet2D *sets, double *rays,
int algorithm)
{
int rval = 0;
log_info("Running benchmark...\n");
double total_initial_time = get_user_time();
stats_printf("cpu_time:\n");
for (int j = 0; j < n_sets; j++)
{
log_debug("Set %d...\n", j);
double set_initial_time = get_user_time();
struct LFreeSet2D *set = &sets[j];
double *pre_m = &PRE_M[j * 4];
double center = CENTER[j];
rval = LFREE_2D_print_set(set);
abort_if(rval, "LFREE_2D_print_set failed");
int wrong_answer = 0;
rval = benchmark_set(algorithm, j, set, rays, pre_m, center, &wrong_answer);
abort_if(rval, "benchmark_set failed");
double set_duration = get_user_time() - set_initial_time;
double avg = (set_duration / N_SAMPLES_PER_SET) * 1000;
if(wrong_answer) avg = 1000000;
stats_printf(" %d: %.8lf\n", j, avg);
log_info(" %3d: %12.3lf ms\n", j, avg);
}
double total_duration = get_user_time() - total_initial_time;
log_info(" %.3lf ms per set \n",
total_duration / (n_sets * N_SAMPLES_PER_SET) * 1000);
if(algorithm == ALGORITHM_MIP)
{
log_info(" %.3lf s spent on LP_create\n", MIP_TIME_CREATE);
log_info(" %.3lf s spent on LP_optimize\n", MIP_TIME_OPTIMIZE);
}
CLEANUP:
return rval;
}
int main(int argc, char **argv)
{
int rval = 0;
double *rays = 0;
struct LFreeSet2D sets[MAX_N_SETS];
rval = parse_args(argc, argv);
if (rval) return 1;
print_header(argc, argv);
if (LOG_FILENAME[0])
{
LOG_FILE = fopen(LOG_FILENAME, "w");
abort_if(!LOG_FILE, "could not open log file");
log_info("Writing log to file: %s\n", LOG_FILENAME);
}
if (STATS_FILENAME[0])
{
STATS_FILE = fopen(STATS_FILENAME, "w");
abort_if(!STATS_FILE, "could not open stats file");
log_info("Writing stats to file: %s\n", STATS_FILENAME);
}
if (WRITE_ANSWERS)
{
N_SAMPLES_PER_SET = 1;
ANSWERS_FILE = fopen(ANSWERS_FILENAME, "w");
abort_if(!ANSWERS_FILE, "could not open answers file");
log_info("Writing answers to file: %s\n", ANSWERS_FILENAME);
}
if (CHECK_ANSWERS)
{
ANSWERS_FILE = fopen(ANSWERS_FILENAME, "r");
abort_if(!ANSWERS_FILE, "could not open answers file");
log_info("Reading answers from file: %s\n", ANSWERS_FILENAME);
}
if(SEED == 0)
{
struct timeval t1;
gettimeofday(&t1, NULL);
SEED = (unsigned int) (t1.tv_usec * t1.tv_sec) % 10000;
}
log_info("Random seed: %u\n", SEED);
srand(SEED);
log_info("Generating %d random rays...\n", N_RAYS);
rays = (double*) malloc(2 * N_RAYS * sizeof(double));
abort_if(!rays, "could not allocate rays");
for (int i = 0; i < N_RAYS; i++)
{
double *ray = &rays[2 * i];
ray[0] = DOUBLE_random(0.0, 1.0);
ray[1] = DOUBLE_random(0.0, 1.0);
}
int algorithm = ALGORITHM_BOUND;
if(SELECT_MIP_ALGORITHM) algorithm = ALGORITHM_MIP;
else if(SELECT_NAIVE_ALGORITHM) algorithm = ALGORITHM_NAIVE;
if(algorithm == ALGORITHM_NAIVE)
{
log_info("Enabling naive algorithm\n");
if(USE_FIXED_BOUNDS)
log_info("Using fixed big M: %d\n", NAIVE_BIG_M);
else
log_info("Enabling bounding boxes\n");
}
else
{
log_info("Enabling bound algorithm\n");
}
log_info("Setting %d samples per set\n", N_SAMPLES_PER_SET);
if(ENABLE_PREPROCESSING)
log_info("Enabling pre-processing\n");
log_info("Reading sets from file...\n");
FILE *sets_file = fopen(SETS_FILENAME, "r");
abort_iff(!sets_file, "could not read file %s", SETS_FILENAME);
int line = 0;
int n_sets = 0;
while(!feof(sets_file))
{
line++;
struct LFreeSet2D *set = &sets[n_sets];
LFREE_2D_init(set, 4, 4, 4);
rval = LFREE_2D_read_next(sets_file, set);
abort_iff(rval, "LFREE_2D_read_next failed (line %d)", line);
if(ENABLE_SHEAR)
{
double m[4] = { 51.0, 5.0, 10.0, 1.0 };
rval = LFREE_2D_transform_set(set, m);
abort_iff(rval, "LFREE_2D_transform_set failed (line %d)", line);
}
double dx = -floor(set->f[0]);
double dy = -floor(set->f[1]);
rval = LFREE_2D_translate_set(set, dx, dy);
abort_iff(rval, "LFREE_2D_translate_set failed (line %d)", line);
if(ENABLE_PREPROCESSING)
{
double *pre_m = &PRE_M[n_sets * 4];
double *center = &CENTER[n_sets];
rval = LFREE_2D_preprocess(set, pre_m, center);
abort_iff(rval, "LFREE_2D_preprocess failed (line %d)", line);
}
rval = LFREE_2D_compute_halfspaces(set);
abort_iff(rval, "LFREE_2D_compute_halfspaces failed (line %d)", line);
rval = LIFTING_2D_verify(set);
if(rval)
{
log_warn(" skipping invalid set on line %d\n", line);
continue;
}
n_sets++;
if(n_sets >= MAX_N_SETS) break;
}
fclose(sets_file);
log_info("Successfully read %d sets\n", n_sets);
rval = benchmark(n_sets, sets, rays, algorithm);
abort_if(rval, "benchmark failed");
log_info("Done.\n");
CLEANUP:
if (LOG_FILE) fclose(LOG_FILE);
if (STATS_FILE) fclose(STATS_FILE);
if (ANSWERS_FILE) fclose(ANSWERS_FILE);
if (rays) free(rays);
return rval;
}

17
lifting/library/CMakeLists.txt
Unescape View File

@ -0,0 +1,17 @@
include_directories(include)
set(COMMON_SOURCES
src/lifting.c
src/lifting-mip.c
include/lifting/lifting-mip.h
include/lifting/lifting.h)
set(TEST_SOURCES
tests/lifting-test.cpp)
add_library(lifting_static ${COMMON_SOURCES})
set_target_properties(lifting_static PROPERTIES OUTPUT_NAME lifting)
target_include_directories (lifting_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
add_executable(lifting-test.run ${COMMON_SOURCES} ${TEST_SOURCES})
target_link_libraries(lifting-test.run gtest_main multirow_static)

32
lifting/library/include/lifting/lifting-mip.h
Unescape View File

@ -0,0 +1,32 @@
/* Copyright (c) 2016 Laurent Poirrier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef LIFTING_MIP_H
#define LIFTING_MIP_H
extern double MIP_TIME_OPTIMIZE;
extern double MIP_TIME_CREATE;
int LIFTING_2D_mip_init();
void LIFTING_2D_mip_cleanup();
int LIFTING_2D_mip(int n_halfspaces,
const double *halfspaces,
const double *ray,
double *value);
#endif

53
lifting/library/include/lifting/lifting.h
Unescape View File

@ -0,0 +1,53 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef LIFTING_H
#define LIFTING_H
#include <multirow/lfree2d.h>
int LIFTING_2D_psi(int n_halfspaces,
const double *halfspaces,
const double *ray,
double *value);
int LIFTING_2D_optimize_continuous(int n_halfspaces,
const double *halfspaces,
double alpha2,
double *alpha1,
double *value);
int LIFTING_2D_lift_fixed(int n_halfspaces,
const double *halfspaces,
const double *ray,
double k1,
double *opt);
int LIFTING_2D_naive(int n_halfspaces,
const double *halfspaces,
const double *ray,
const int *lb,
const int *ub,
double *value);
int LIFTING_2D_bound(int n_halfspaces,
const double *halfspaces,
const double *ray,
double *value);
int LIFTING_2D_verify(struct LFreeSet2D *set);
#endif //LIFTING_H

162
lifting/library/src/lifting-mip.c
Unescape View File

@ -0,0 +1,162 @@
/* Copyright (c) 2016 Laurent Poirrier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <multirow/util.h>
#include <multirow/lp.h>
#include <lifting/lifting-mip.h>
double MIP_TIME_OPTIMIZE = 0;
double MIP_TIME_CREATE = 0;
/**
* Returns non-zero if the cplex_env part of struct LP is initialized.
*
* @param lp pointer to a struct LP
*/
int LP_is_open(struct LP *lp)
{
return (lp->cplex_env != 0);
}
/**
* Frees the cplex_lp part of struct LP, leaving the environment open.
*
* @param lp pointer to a struct LP
*/
void LP_destroy(struct LP *lp)
{
if(!lp) return;
if(!lp->cplex_env) return;
if(lp->cplex_lp)
CPXfreeprob(lp->cplex_env, &(lp->cplex_lp));
lp->cplex_lp = 0;
}
/**
* Static data for LIFTING_2D_mip().
*
* Warning: Because of this, LIFTING_2D_mip() is not thread safe.
*/
static struct LP LIFTING_2D_mip_lp = {0, 0};
/**
* Opens the cplex environment of the struct LP used by LIFTING_2D_mip().
*
* This is called automatically by LIFTING_2D_mip().
*/
int LIFTING_2D_mip_init()
{
struct LP *lp = &LIFTING_2D_mip_lp;
if(!LP_is_open(lp))
{
if(LP_open(lp)) return (-1);
}
return (0);
}
/**
* Closes the cplex environment of the struct LP used by LIFTING_2D_mip().
*
* This function should be called by the user of LIFTING_2D_mip(), if
* cleaning up the environment is desired.
*/
void LIFTING_2D_mip_cleanup()
{
struct LP *lp = &LIFTING_2D_mip_lp;
LP_free(lp);
}
/*
* Computes the lifting coefficient of a ray by formulating the lifting
* problem as a MIP.
*
* @param[in] n_halfspaces number of facets of the lattice-free set used
* @param[in] halfspaces description of the lattice-free set used
* @param[in] ray ray to lift
* @param[out] value lifting coefficient
*/
int LIFTING_2D_mip(int n_halfspaces,
const double *halfspaces,
const double *ray,
double *value)
{
struct LP *lp = &LIFTING_2D_mip_lp;
int rval = 0;
double initial_time = get_user_time();
rval = LIFTING_2D_mip_init();
abort_if(rval, "LIFTING_2D_mip_init failed");
rval = LP_create(lp, "lifting2d");
abort_if(rval, "LP_create failed");
rval = LP_new_col(lp, 1.0, 0.0, MILP_INFINITY, 'C');
rval |= LP_new_col(lp, 0.0, -MILP_INFINITY, MILP_INFINITY, 'I');
rval |= LP_new_col(lp, 0.0, -MILP_INFINITY, MILP_INFINITY, 'I');
abort_if(rval, "LP_new_col failed");
double rhs;
char sense = 'G';
int beg = 0;
int ind[3] = {0, 1, 2};
double val[3];
val[0] = 1.0;
for(int i = 0; i < n_halfspaces; i++)
{
double a0 = halfspaces[i * 2 + 0];
double a1 = halfspaces[i * 2 + 1];
rhs = a0 * ray[0] + a1 * ray[1];
val[1] = -a0;
val[2] = -a1;
rval |= LP_add_rows(lp, 1, 3, &rhs, &sense, &beg, ind, val);
}
MIP_TIME_CREATE += get_user_time() - initial_time;
abort_if(rval, "LP_add_rows failed");
int infeasible;
initial_time = get_user_time();
rval = LP_optimize(lp, &infeasible);
abort_if(rval, "LP_optimize failed");
abort_if(infeasible, "LIFTING_2D_mip infeasible");
MIP_TIME_OPTIMIZE += get_user_time() - initial_time;
double obj;
rval = LP_get_obj_val(lp, &obj);
abort_if(rval, "LP_get_obj_val failed");
*value = obj;
CLEANUP:
LP_destroy(lp);
return (rval);
}

286
lifting/library/src/lifting.c
Unescape View File

@ -0,0 +1,286 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define _XOPEN_SOURCE 700
#include <math.h>
#include <stdio.h>
#include <multirow/double.h>
#include <multirow/util.h>
#include <multirow/lp.h>
#include <multirow/lfree2d.h>
#include <lifting/lifting.h>
/**
* Verifies if the set is well formed.
*/
int LIFTING_2D_verify(struct LFreeSet2D *set)
{
int rval = 0;
abort_if(set->n_halfspaces == 0, "Halfspaces not found");
for(int i = 0; i < set->n_lattice_points; i++)
{
double *t = &set->lattice_points[2 * i];
double r[2] = {t[0] - set->f[0],
t[1] - set->f[1]};
double value;
rval = LIFTING_2D_psi(set->n_halfspaces, set->halfspaces, r, &value);
abort_if(rval, "LIFTING_2D_psi failed");
double delta = fabs(value - 1);
if(delta > 0.0001)
{
log_debug("Lattice point (%.2lf, %.2lf) is "
"not on the boundary (delta=%.6lf)", t[0], t[1], delta);
rval = 1;
goto CLEANUP;
}
}
CLEANUP:
return rval;
}
int LIFTING_2D_psi(int n_halfspaces,
const double *halfspaces,
const double *ray,
double *value)
{
int rval = 0;
*value = -INFINITY;
for(int i = 0; i < n_halfspaces; i++)
{
const double *h = &halfspaces[2 * i];
double v = ray[0] * h[0] + ray[1] * h[1];
*value = DOUBLE_max(v, *value);
}
CLEANUP:
return rval;
}
int LIFTING_2D_optimize_continuous(int n_halfspaces,
const double *halfspaces,
double alpha2,
double *alpha1,
double *value)
{
int rval = 0;
*value = -INFINITY;
for(int r = 0; r < n_halfspaces; r++)
{
for(int s = r+1; s < n_halfspaces; s++)
{
const double *ar = &halfspaces[r * 2];
const double *as = &halfspaces[s * 2];
double xr = ar[0];
double xs = as[0];
if(DOUBLE_eq(xr, xs)) continue;
double br = - xs / (xr - xs);
double bs = xr / (xr - xs);
if(!DOUBLE_geq(br, 0)) continue;
if(!DOUBLE_geq(bs, 0)) continue;
double yr = ar[1] * alpha2;
double ys = as[1] * alpha2;
double obj = yr * br + ys * bs;
if(DOUBLE_geq(obj, *value))
{
*value = obj;
if(DOUBLE_neq(xr, 0))
*alpha1 = - (yr - obj) / xr;
else
*alpha1 = - (ys - obj) / xs;
}
}
}
CLEANUP:
return rval;
}
int LIFTING_2D_lift_fixed(int n_halfspaces,
const double *halfspaces,
const double *ray,
double k1,
double *opt)
{
int rval = 0;
double k0;
double value;
rval = LIFTING_2D_optimize_continuous(n_halfspaces, halfspaces,
ray[1] + k1, &k0, &value);
abort_if(rval, "LIFTING_2D_optimize_continuous failed");
double delta = k0 - ray[0];
double r_ceil[2] = { ray[0] + ceil(delta), ray[1] + k1 };
double r_floor[2] = { ray[0] + floor(delta), ray[1] + k1 };
log_debug(" delta=%.6lf value=%.6lf\n", delta, value);
log_debug(" r_ceil=%12.6lf %12.6lf\n", r_ceil[0], r_ceil[1]);
log_debug(" r_floor=%12.6lf %12.6lf\n", r_floor[0], r_floor[1]);
double value_ceil, value_floor;
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r_ceil, &value_ceil);
abort_if(rval, "LIFTING_2D_psi failed");
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r_floor, &value_floor);
abort_if(rval, "LIFTING_2D_psi failed");
*opt = min(value_ceil, value_floor);
CLEANUP:
return rval;
}
int LIFTING_2D_naive(int n_halfspaces,
const double *halfspaces,
const double *ray,
const int *lb,
const int *ub,
double *value)
{
int rval = 0;
*value = INFINITY;
int best_k0 = 0;
int best_k1 = 0;
int margin = 10;
for(int k0 = lb[0] - margin; k0 <= ub[0] + margin; k0++)
{
for(int k1 = lb[1] - margin; k1 <= ub[1] + margin; k1++)
{
double q[2] = { ray[0] + k0, ray[1] + k1 };
double value_q;
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, q, &value_q);
abort_if(rval, "LIFTING_2D_ps failed");
if(value_q < *value)
{
best_k0 = k0;
best_k1 = k1;
*value = value_q;
log_debug(" k=%6d %6d value=%12.6lf\n", k0, k1, *value);
}
}
}
// log_debug(" best_k=(%d %d) value=%.6lf\n", best_k0, best_k1, *value);
CLEANUP:
return rval;
}
int LIFTING_2D_bound(int n_halfspaces,
const double *halfspaces,
const double *ray,
double *value)
{
int rval = 0;
double eta_star, eta_plus, eta_minus;
double xi_plus, xi_minus;
double m_plus, m_minus;
double ignored;
int k1 = 1;
rval = LIFTING_2D_lift_fixed(n_halfspaces, halfspaces, ray, 0, &eta_star);
abort_if(rval, "LIFTING_2D_lift_fixed failed");
rval = LIFTING_2D_optimize_continuous(n_halfspaces, halfspaces, 1, &ignored,
&xi_plus);
abort_if(rval, "LIFTING_2D_optimize_continuous failed");
rval = LIFTING_2D_optimize_continuous(n_halfspaces, halfspaces, -1,
&ignored, &xi_minus);
abort_if(rval, "LIFTING_2D_optimize_continuous failed");
m_plus = m_minus = INFINITY;
log_debug("Level 0:\n");
log_debug(" eta star = %.6lf\n", eta_star);
log_debug(" next slack plus = %.6lf\n", k1 + ray[1] - m_plus);
log_debug(" next slack minus = %.6lf\n", k1 - ray[1] - m_minus);
log_debug(" xi plus = %.6lf\n", xi_plus);
log_debug(" xi minus = %.6lf\n", xi_minus);
while((k1 <= fabs(ray[1])) || (k1 + ray[1] <= m_plus) || (k1 - ray[1] <= m_minus))
{
log_debug("Level %d:\n", k1);
double h_plus, h_minus;
rval = LIFTING_2D_optimize_continuous(n_halfspaces, halfspaces, ray[1]
+ k1, &ignored, &h_plus);
abort_if(rval, "LIFTING_2D_optimize_continuous failed");
rval = LIFTING_2D_optimize_continuous(n_halfspaces, halfspaces, ray[1]
- k1, &ignored, &h_minus);
abort_if(rval, "LIFTING_2D_optimize_continuous failed");
rval = LIFTING_2D_lift_fixed(n_halfspaces, halfspaces, ray, k1,
&eta_plus);
abort_if(rval, "LIFTING_2D_lift_fixed failed");
rval = LIFTING_2D_lift_fixed(n_halfspaces, halfspaces, ray, -k1,
&eta_minus);
abort_if(rval, "LIFTING_2D_lift_fixed failed");
eta_star = fmin(eta_star, fmin(eta_plus, eta_minus));
m_plus = ceil(eta_star / xi_plus);
m_minus = ceil(eta_star / xi_minus);
log_debug(" eta plus = %.6lf\n", eta_plus);
log_debug(" eta minus = %.6lf\n", eta_minus);
log_debug(" eta star = %.6lf\n", eta_star);
log_debug(" h minus = %.6lf\n", h_minus);
log_debug(" h plus = %.6lf\n", h_plus);
k1++;
log_debug(" next slack plus = %.6lf\n", k1 + ray[1] - m_plus);
log_debug(" next slack minus = %.6lf\n", k1 - ray[1] - m_minus);
}
log_debug("Done\n");
*value = eta_star;
CLEANUP:
return rval;
}

1
lifting/library/tests/fixtures/quads.txt vendored
Unescape View File

@ -0,0 +1 @@
0.5 0.5 4 -0.5 0.5 0.5 1.5 1.5 0.5 0.5 -0.5 4 0 0 0 1 1 1 1 0

2
lifting/library/tests/fixtures/triangles.txt vendored
Unescape View File

@ -0,0 +1,2 @@
0.7875 1.7875 3 1 2 0.4875 2.0875 -0.0625 -16.0625 3 1 2 0 -14 0 -15
0.5 0.5 3 0 0 2 0 0 2 3 0 1 1 0 1 1

715
lifting/library/tests/lifting-test.cpp
Unescape View File

@ -0,0 +1,715 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <gtest/gtest.h>
extern "C" {
#include <multirow/util.h>
#include <multirow/double.h>
#include <multirow/lfree2d.h>
#include <lifting/lifting.h>
#include <lifting/lifting-mip.h>
#include <math.h>
}
#define E 1e-6
TEST(Lifting2DTest, lift_fixed_test_2)
{
int rval = 0;
double opt;
int n_halfspaces = 3;
double halfspaces[] = {
-1.12796209, -2.25592417,
2.38125329, 11.19606810,
7.35264174, 6.22467965,
};
double r[] = { 0.48231629, 0.70551355 };
rval = LIFTING_2D_lift_fixed(n_halfspaces, halfspaces, r, -1, &opt);
abort_if(rval, "LIFTING_lift_fixed failed");
EXPECT_NEAR(opt, 1.24826670, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, psi_test)
{
int rval = 0;
int n_halfspaces = 5;
double halfspaces[] = {
-1 / 2.0, 0,
-1 / 2.0, -1 / 2.0,
-1 / 4.0, 1 / 2.0,
-1 / 6.0, -5 / 6.0,
5 / 16.0, 1 / 8.0
};
double r1[] = { 0, 0 };
double r2[] = { 2, 3 };
double r3[] = { 1 / 3.0, 1 / 2.0 };
double r4[] = { -4.0, 0 };
double r5[] = { 1, -1/2.0 };
double r6[] = { -1/2.0, 1/2.0 };
double r7[] = { 1/2.0, 1/2.0 };
double value;
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r1, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 0.0, E);
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r2, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 1.0, E);
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r3, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 1 / 6.0, E);
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r4, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 2.0, E);
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r5, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 1 / 4.0, E);
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r6, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 0.375, E);
rval = LIFTING_2D_psi(n_halfspaces, halfspaces, r7, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 7 / 32.0, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, psi_test_2)
{
int rval = 0;
LFreeSet2D set;
set.f[0] = 0;
set.f[1] = 0.359;
double vertices[] = {
-1.700870, -0.796700,
0.068700, 0.032200,
1.106900, 1.111100,
-0.131990, 0.986800,
};
double lattice_points[] = {
0.000000, 0.000000,
1.000000, 1.000000,
0.000000, 1.000000,
-1.000000, 0.000000,
};
int ub[] = { 10, 10 };
int lb[] = { -12, -12 };
double r0[] = { 0.37, 0.19 };
double value;
rval = LFREE_2D_init(&set, 100, 100, 100);
abort_if(rval, "LFREE_2D_init failed");
set.n_vertices = 4;
set.vertices = vertices;
set.n_lattice_points = 4;
set.lattice_points = lattice_points;
rval = LFREE_2D_compute_halfspaces(&set);
abort_if(rval, "LFREE_2D_compute_halfspaces failed");
rval = LFREE_2D_print_set(&set);
abort_if(rval, "LFREE_2D_print_set failed");
rval = LIFTING_2D_naive(set.n_halfspaces, set.halfspaces, r0, lb, ub, &value);
abort_if(rval, "LIFTING_2D_psi failed");
EXPECT_NEAR(value, 0.48846868, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, naive_test_1)
{
int rval = 0;
int n_halfspaces = 5;
double halfspaces[] = {
-1 / 2.0, 0,
-1 / 2.0, -1 / 2.0,
-1 / 4.0, 1 / 2.0,
-1 / 6.0, -5 / 6.0,
5 / 16.0, 1 / 8.0
};
double r1[] = { 1/2.0, 1/2.0 };
double r2[] = { 0, 1/2.0 };
double r3[] = { 1/3.0, 2/3.0 };
double r4[] = { 2/5.0, 3/7.0 };
int lb[] = { -10, -10 };
int ub[] = { 10, 10 };
double value;
rval = LIFTING_2D_naive(n_halfspaces, halfspaces, r1, lb, ub, &value);
abort_if(rval, "LIFTING_2D_naive failed");
EXPECT_NEAR(value, 0.21875, E);
rval = LIFTING_2D_naive(n_halfspaces, halfspaces, r2, lb, ub, &value);
abort_if(rval, "LIFTING_2D_naive failed");
EXPECT_NEAR(value, 0.25, E);
rval = LIFTING_2D_naive(n_halfspaces, halfspaces, r3, lb, ub, &value);
abort_if(rval, "LIFTING_2D_naive failed");
EXPECT_NEAR(value, 0.222222, E);
rval = LIFTING_2D_naive(n_halfspaces, halfspaces, r4, lb, ub, &value);
abort_if(rval, "LIFTING_2D_naive failed");
EXPECT_NEAR(value, 0.178571, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, naive_test_2)
{
int rval = 0;
int n_halfspaces = 3;
double halfspaces[] = {
-1.12796209, -2.25592417,
2.38125329, 11.19606810,
7.35264174, 6.22467965,
};
double r[] = { 0.48231629, 0.70551355 };
int lb[] = { -50, -50 };
int ub[] = { 50, 50 };
double value;
rval = LIFTING_2D_naive(n_halfspaces, halfspaces, r, lb, ub, &value);
abort_if(rval, "LIFTING_2D_naive failed");
EXPECT_NEAR(value, 1.24826671, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, bound_test_1)
{
int rval = 0;
int n_halfspaces = 5;
double halfspaces[] = {
-1 / 2.0, 0,
-1 / 2.0, -1 / 2.0,
-1 / 4.0, 1 / 2.0,
-1 / 6.0, -5 / 6.0,
5 / 16.0, 1 / 8.0
};
double r1[] = { 1/2.0, 1/2.0 };
double r2[] = { 0, 1/2.0 };
double r3[] = { 1/3.0, 2/3.0 };
double r4[] = { 2/5.0, 3/7.0 };
double value;
rval = LIFTING_2D_bound(n_halfspaces, halfspaces, r1, &value);
abort_if(rval, "LIFTING_2D_bound failed");
EXPECT_NEAR(value, 0.21875, E);
rval = LIFTING_2D_bound(n_halfspaces, halfspaces, r2, &value);
abort_if(rval, "LIFTING_2D_bound failed");
EXPECT_NEAR(value, 0.25, E);
rval = LIFTING_2D_bound(n_halfspaces, halfspaces, r3, &value);
abort_if(rval, "LIFTING_2D_bound failed");
EXPECT_NEAR(value, 0.222222, E);
rval = LIFTING_2D_bound(n_halfspaces, halfspaces, r4, &value);
abort_if(rval, "LIFTING_2D_bound failed");
EXPECT_NEAR(value, 0.178571, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, bound_test_2)
{
int rval = 0;
int n_halfspaces = 3;
double halfspaces[] = {
-1.12796209, -2.25592417,
2.38125329, 11.19606810,
7.35264174, 6.22467965,
};
double r[] = { 0.48231629, 0.70551355 };
double value;
rval = LIFTING_2D_bound(n_halfspaces, halfspaces, r, &value);
abort_if(rval, "LIFTING_2D_bound failed");
EXPECT_NEAR(value, 1.24826671, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, mip_test_1)
{
int rval = 0;
int n_halfspaces = 5;
double halfspaces[] = {
-1 / 2.0, 0,
-1 / 2.0, -1 / 2.0,
-1 / 4.0, 1 / 2.0,
-1 / 6.0, -5 / 6.0,
5 / 16.0, 1 / 8.0
};
double r1[] = { 1/2.0, 1/2.0 };
double r2[] = { 0, 1/2.0 };
double r3[] = { 1/3.0, 2/3.0 };
double r4[] = { 2/5.0, 3/7.0 };
double value;
rval = LIFTING_2D_mip(n_halfspaces, halfspaces, r1, &value);
abort_if(rval, "LIFTING_2D_mip failed");
EXPECT_NEAR(value, 0.21875, E);
rval = LIFTING_2D_mip(n_halfspaces, halfspaces, r2, &value);
abort_if(rval, "LIFTING_2D_mip failed");
EXPECT_NEAR(value, 0.25, E);
rval = LIFTING_2D_mip(n_halfspaces, halfspaces, r3, &value);
abort_if(rval, "LIFTING_2D_mip failed");
EXPECT_NEAR(value, 0.222222, E);
rval = LIFTING_2D_mip(n_halfspaces, halfspaces, r4, &value);
abort_if(rval, "LIFTING_2D_mip failed");
EXPECT_NEAR(value, 0.178571, E);
CLEANUP:
if(rval) FAIL();
}
TEST(Lifting2DTest, verify_test_1)
{
int rval = 0;
double vertices[] = {
0.361588, -0.411851,
2.550110, 1.000000,
1.009740, 1.000000,
-0.836467, 0.952742,
};
double lattice_points[] = {
1.000000, 0.000000,
2.000000, 1.000000,
1.000000, 1.000000,
0.000000, 0.000000,
};
LFreeSet2D set;
rval = LFREE_2D_init(&set, 4, 4, 4);
abort_if(rval, "LFREE_2D_init failed");
set.n_vertices = 4;
set.vertices = vertices;
set.n_lattice_points = 4;
set.lattice_points = lattice_points;
set.f[0] = 0.918857;
set.f[1] = 0.952742;
rval = LFREE_2D_compute_halfspaces(&set);
abort_if(rval, "LFREE_2D_compute_halfspaces failed");
rval = LIFTING_2D_verify(&set);
EXPECT_EQ(rval, 1);
rval = 0;
CLEANUP:
if(rval) FAIL();
}
TEST(LFreeSetTest, read_next_test)
{
int rval = 0;
FILE *triangles = fopen("../lifting/library/tests/fixtures/triangles.txt", "r");
abort_if(!triangles, "could not read triangles.txt");
LFreeSet2D set;
LFREE_2D_init(&set, 100, 100, 100);
rval = LFREE_2D_read_next(triangles, &set);
abort_if(rval, "LFREE_2D_read_next failed");
EXPECT_NEAR(set.f[0], 0.7875, E);
EXPECT_NEAR(set.f[1], 1.7875, E);
EXPECT_NEAR(set.vertices[0], 1, E);
EXPECT_NEAR(set.vertices[1], 2, E);
EXPECT_NEAR(set.vertices[2], 0.4875, E);
EXPECT_NEAR(set.vertices[3], 2.0875, E);
EXPECT_NEAR(set.vertices[4], -0.0625, E);
EXPECT_NEAR(set.vertices[5], -16.0625, E);
EXPECT_NEAR(set.lattice_points[0], 1, E);
EXPECT_NEAR(set.lattice_points[1], 2, E);
EXPECT_NEAR(set.lattice_points[2], 0, E);
EXPECT_NEAR(set.lattice_points[3], -14, E);
EXPECT_NEAR(set.lattice_points[4], 0, E);
EXPECT_NEAR(set.lattice_points[5], -15, E);
rval = LFREE_2D_compute_halfspaces(&set);
abort_if(rval, "LFREE_2D_compute_halfspaces failed");
EXPECT_NEAR(set.halfspaces[0], 0.686274, E);
EXPECT_NEAR(set.halfspaces[1], 4.019607, E);
EXPECT_NEAR(set.halfspaces[2], -3.235294, E);
EXPECT_NEAR(set.halfspaces[3], 0.098039, E);
EXPECT_NEAR(set.halfspaces[4], 5.0, E);
EXPECT_NEAR(set.halfspaces[5], -0.294117, E);
CLEANUP:
if(rval) FAIL();
}
TEST(LFreeSetTest, read_next_quadrilateral_test)
{
int rval = 0;
FILE *file = fopen("../lifting/library/tests/fixtures/quads.txt", "r");
abort_if(!file, "could not read quads.txt");
LFreeSet2D set;
LFREE_2D_init(&set, 100, 100, 100);
rval = LFREE_2D_read_next(file, &set);
abort_if(rval, "LFREE_2D_read_next failed");
EXPECT_NEAR(set.f[0], 0.5, E);
EXPECT_NEAR(set.f[1], 0.5, E);
EXPECT_NEAR(set.vertices[0], -0.5, E);
EXPECT_NEAR(set.vertices[1], 0.5, E);
EXPECT_NEAR(set.vertices[2], 0.5, E);
EXPECT_NEAR(set.vertices[3], 1.5, E);
EXPECT_NEAR(set.vertices[4], 1.5, E);
EXPECT_NEAR(set.vertices[5], 0.5, E);
EXPECT_NEAR(set.vertices[6], 0.5, E);
EXPECT_NEAR(set.vertices[7], -0.5, E);
EXPECT_NEAR(set.lattice_points[0], 0, E);
EXPECT_NEAR(set.lattice_points[1], 0, E);
EXPECT_NEAR(set.lattice_points[2], 0, E);
EXPECT_NEAR(set.lattice_points[3], 1, E);
EXPECT_NEAR(set.lattice_points[4], 1, E);
EXPECT_NEAR(set.lattice_points[5], 1, E);
EXPECT_NEAR(set.lattice_points[6], 1, E);
EXPECT_NEAR(set.lattice_points[7], 0, E);
rval = LFREE_2D_compute_halfspaces(&set);
abort_if(rval, "LFREE_2D_compute_halfspaces failed");
EXPECT_NEAR(set.halfspaces[0], -1, E);
EXPECT_NEAR(set.halfspaces[1], 1, E);
EXPECT_NEAR(set.halfspaces[2], 1, E);
EXPECT_NEAR(set.halfspaces[3], 1, E);
EXPECT_NEAR(set.halfspaces[4], 1, E);
EXPECT_NEAR(set.halfspaces[5], -1, E);
EXPECT_NEAR(set.halfspaces[6], -1, E);
EXPECT_NEAR(set.halfspaces[7], -1, E);
CLEANUP:
if(rval) FAIL();
}
//TEST(LFreeSetTest, get_bounding_box_test)
//{
// int rval = 0;
//
// FILE *triangles = fopen("../tests/fixtures/triangles.txt", "r");
// abort_if(!triangles, "could not read triangles.txt");
//
// LFreeSet2D set;
// LFREE_2D_init(&set, 100, 100, 100);
//
// int ub[2], lb[2];
//
// rval = LFREE_2D_read_next(triangles, &set);
// abort_if(rval, "LFREE_2D_read_next failed");
//
// rval = LFREE_2D_get_bounding_box(&set, lb, ub);
// abort_if(rval, "LFREE_2D_get_bounding_box failed");
//
// EXPECT_NEAR(lb[0], -1.0, E);
// EXPECT_NEAR(ub[0], 1.0, E);
//
// EXPECT_NEAR(lb[1], -18.0, E);
// EXPECT_NEAR(ub[1], 2, E);
//
//CLEANUP:
// if(rval) FAIL();
//}
TEST(LFreeSetTest, preprocess_test)
{
int rval = 0;
double vertices[] = {
100/11.0 ,167/33.0,
10/11.0 , 97/33.0,
20/11.0 , 39/11.0
};
double lattice_points[] = {
1, 3,
4, 4,
5, 4
};
double pre_m[4];
double center;
struct LFreeSet2D set;
set.f[0] = 9 / 2.0;
set.f[1] = 4.0;
set.n_vertices = 3;
set.vertices = vertices;
set.n_lattice_points = 3;
set.lattice_points = lattice_points;
rval = LFREE_2D_preprocess(&set, pre_m, &center);
abort_if(rval, "LFREE_2D_preprocess failed");
EXPECT_NEAR(set.f[0], 1 / 2.0, E);
EXPECT_NEAR(set.f[1], 1 / 2.0, E);
EXPECT_NEAR(set.vertices[0], 21 / 11.0, E);
EXPECT_NEAR(set.vertices[1], 5 / 33.0, E);
EXPECT_NEAR(set.vertices[2], 1 / 11.0, E);
EXPECT_NEAR(set.vertices[3], -5 / 33.0, E);
EXPECT_NEAR(set.vertices[4], -9 / 11.0, E);
EXPECT_NEAR(set.vertices[5], 15 / 11.0, E);
EXPECT_NEAR(center, 20 / 33.0, E);
CLEANUP:
if(rval) FAIL();
}TEST(LFreeSetTest, preprocess_test_3)
{
int rval = 0;
double vertices[] = {
1.0, -4.0,
1.0, 2.0,
-0.2, 0.8
};
double lattice_points[] = {
1, 0,
0, 1,
0, 0
};
double pre_m[4];
double center;
struct LFreeSet2D set;
set.f[0] = 1 / 2.0;
set.f[1] = 1 / 2.0;
set.n_vertices = 3;
set.vertices = vertices;
set.n_lattice_points = 3;
set.lattice_points = lattice_points;
rval = LFREE_2D_preprocess(&set, pre_m, &center);
abort_if(rval, "LFREE_2D_preprocess failed");
EXPECT_NEAR(set.vertices[0], -4, E);
EXPECT_NEAR(set.vertices[1], 1, E);
EXPECT_NEAR(set.vertices[2], 2, E);
EXPECT_NEAR(set.vertices[3], 1, E);
EXPECT_NEAR(set.vertices[4], 0.8, E);
EXPECT_NEAR(set.vertices[5], -0.2, E);
EXPECT_NEAR(center, 0.4, E);
CLEANUP:
if(rval) FAIL();
}
TEST(LFreeSetTest, preprocess_test_2)
{
int rval = 0;
double vertices[] = {
22, 69/7.0,
-3, -11/7.0,
-8, -26/7.0
};
double lattice_points[] = {
-4, -2,
-2, -1,
7, 3
};
double r[] = { 2/3.0, 2/6.0 };
double r0, r1;
double value;
double pre_m[4];
struct LFreeSet2D set;
set.f[0] = 2 / 3.0;
set.f[1] = 1 / 6.0;
rval = LFREE_2D_init(&set, 10, 10, 10);
abort_if(rval, "LFREE_2D_init failed");
set.n_vertices = 3;
set.vertices = vertices;
set.n_lattice_points = 3;
set.lattice_points = lattice_points;
rval = LFREE_2D_print_set(&set);
abort_if(rval, "LFREE_2D_print_set failed");
rval = LFREE_2D_compute_halfspaces(&set);
abort_if(rval, "LFREE_2D_compute_halfspaces failed");
rval = LIFTING_2D_bound(set.n_halfspaces, set.halfspaces, r, &value);
abort_if(rval, "LIFTING_2D_bound failed");
log_debug("%.6lf\n", value);
CLEANUP:
if(rval) FAIL();
}
TEST(LFreeSetTest, preprocess_test_4)
{
int rval = 0;
double vertices[] = {
-5/2.0, -1,
17/4.0, 2,
7/2.0, 2,
-13/4.0, -1
};
double lattice_points[] = {
2, 1,
-3, -1,
-1, 0,
4, 2
};
double pre_m[4];
double center;
struct LFreeSet2D set;
set.f[0] = 1 / 2.0;
set.f[1] = 1 / 2.0;
set.n_vertices = 4;
set.vertices = vertices;
set.n_lattice_points = 4;
set.lattice_points = lattice_points;
rval = LFREE_2D_preprocess(&set, pre_m, &center);
abort_if(rval, "LFREE_2D_preprocess failed");
LFREE_2D_print_set(&set);
EXPECT_NEAR(set.vertices[0], -1, E);
EXPECT_NEAR(set.vertices[1], 1/2.0, E);
EXPECT_NEAR(set.vertices[2], 1/2.0, E);
EXPECT_NEAR(set.vertices[3], -1/4.0, E);
EXPECT_NEAR(set.vertices[4], 2, E);
EXPECT_NEAR(set.vertices[5], 1/2.0, E);
EXPECT_NEAR(set.vertices[6], 1/2.0, E);
EXPECT_NEAR(set.vertices[7], 5/4.0, E);
EXPECT_NEAR(pre_m[0], -2, E);
EXPECT_NEAR(pre_m[1], -1, E);
EXPECT_NEAR(pre_m[2], 5, E);
EXPECT_NEAR(pre_m[3], 2, E);
EXPECT_NEAR(center, 0.5, E);
CLEANUP:
if(rval) FAIL();
}

32
multirow/CMakeLists.txt
Unescape View File

@ -0,0 +1,32 @@
set(COMMON_SOURCES
src/cg.c
src/double.c
src/geometry.c
src/lfree2d.c
src/lp.c
src/mir.c
src/util.c
src/rational.c
src/stats.c
include/multirow/cg.h
include/multirow/double.h
include/multirow/geometry.h
include/multirow/lfree2d.h
include/multirow/lp.h
include/multirow/rational.h
include/multirow/stats.h
include/multirow/params.h
include/multirow/util.h)
set(TEST_SOURCES
tests/double-test.cpp
tests/cg-test.cpp
tests/geometry-test.cpp)
add_library(multirow_static ${COMMON_SOURCES})
set_target_properties(multirow_static PROPERTIES OUTPUT_NAME lifting)
target_link_libraries(multirow_static ${CPLEX_LIBRARIES})
target_include_directories (multirow_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
add_executable(multirow-test.run ${COMMON_SOURCES} ${TEST_SOURCES})
target_link_libraries(multirow-test.run ${CPLEX_LIBRARIES} gtest_main)

95
multirow/include/multirow/cg.h
Unescape View File

@ -0,0 +1,95 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_CG_H
#define MULTIROW_CG_H
#include "lp.h"
struct CG
{
struct LP *lp;
struct Row **tableau_rows;
int *cstat;
int *rstat;
double *ub;
double *lb;
int nrows;
int ncols;
char *column_types;
double last_obj_value;
double *integral_solution;
double *basic_solution;
double *current_solution;
};
typedef int (*SingleRowGeneratorCallback)(const struct Row *row,
char *column_types,
struct Row *cut);
typedef int (*MultirowGeneratorCallback)(int nrows,
struct Row **rows,
char *column_types,
struct Row *cut);
int CG_init(struct LP *lp,
char *column_types,
struct CG *cg);
void CG_free(struct CG *cg);
int CG_add_single_row_cuts(struct CG *cg,
SingleRowGeneratorCallback generate);
int CG_add_multirow_cuts(struct CG *cg,
int nrows,
MultirowGeneratorCallback generate);
int CG_set_integral_solution(struct CG *cg,
double *valid_solution);
int CG_set_basic_solution(struct CG *cg,
double *basic_solution);
int CG_extract_rays_from_rows(int nrows,
struct Row **rows,
double *rays,
int *nrays,
int *variable_to_ray,
double *ray_scale,
int *indices,
int *nz);
int CG_boost_variable(int var,
double factor,
int nrows,
double *rays,
int *variable_to_ray,
double *ray_scale,
int *indices,
int nz);
int CG_find_ray(int dim,
const double *rays,
int nrays,
const double *r,
int *found,
double *scale,
int *index);
#endif //MULTIROW_CG_H

50
multirow/include/multirow/double.h
Unescape View File

@ -0,0 +1,50 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_DOUBLE_H
#define MULTIROW_DOUBLE_H
#include <multirow/rational.h>
int DOUBLE_sgn(double a);
int DOUBLE_cmp(double a,
double b);
int DOUBLE_geq(double a,
double b);
int DOUBLE_leq(double a,
double b);
int DOUBLE_eq(double a,
double b);
int DOUBLE_neq(double a,
double b);
double DOUBLE_max(double a,
double b);
int DOUBLE_to_rational(double a,
long max_den,
Rational r);
#define DOUBLE_iszero(a) (fabs(a) < EPSILON)
double DOUBLE_random(double min, double max);
#endif //MULTIROW_DOUBLE_H

30
multirow/include/multirow/geometry.h
Unescape View File

@ -0,0 +1,30 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_GEOMETRY_H
#define MULTIROW_GEOMETRY_H
int scale_vector_to_line(const double *ray1,
const double *ray2,
const double *p,
double *lambda);
int chull_2d(const double *original_points,
int npoints,
double *vertices,
int *nvertices);
#endif //MULTIROW_GEOMETRY_H

58
multirow/include/multirow/lfree2d.h
Unescape View File

@ -0,0 +1,58 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef LFREE_2D_H
#define LFREE_2D_H
struct LFreeSet2D
{
double f[2];
int n_vertices;
double *vertices;
int n_lattice_points;
double *lattice_points;
int n_halfspaces;
double *halfspaces;
};
int LFREE_2D_init(struct LFreeSet2D *set,
int n_vertices,
int n_lattice_points,
int max_n_halfspaces);
void LFREE_2D_free(struct LFreeSet2D *set);
int LFREE_2D_read_next(FILE *file, struct LFreeSet2D *set);
int LFREE_2D_compute_halfspaces(struct LFreeSet2D *set);
int LFREE_2D_preprocess(struct LFreeSet2D *set, double *m, double *center);
int LFREE_2D_transform_set(struct LFreeSet2D *set, const double *m);
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);
#endif //LFREE_2D_H

157
multirow/include/multirow/lp.h
Unescape View File

@ -0,0 +1,157 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _PROJECT_LP_H_
#define _PROJECT_LP_H_
#include <ilcplex/cplex.h>
#include "params.h"
#define MILP_INTEGER 0
#define MILP_CONTINUOUS 1
#define MILP_INFINITY CPX_INFBOUND
struct LP
{
CPXENVptr cplex_env;
CPXLPptr cplex_lp;
};
struct Row
{
int nz;
int head;
double pi_zero;
double *pi;
int *indices;
};
static const int MAX_NAME_LENGTH = 100;
int LP_open(struct LP *lp);
int LP_create(struct LP *lp,
const char *name);
int LP_write(struct LP *lp,
const char *fname);
void LP_free(struct LP *lp);
void LP_free_row(struct Row *row);
int LP_new_row(struct LP *lp,
char sense,
double rhs);
int LP_new_col(struct LP *lp,
double obj,
double lb,
double ub,
char type);
int LP_add_row(struct LP *lp,
struct Row *row);
int LP_add_rows(struct LP *lp,
int newrows,
int newnz,
double *rhs,
char *sense,
int *rmatbeg,
int *rmatind,
double *rmatval);
int LP_add_cols(struct LP *lp,
int newcols,
int newnz,
double *obj,
int *cmatbeg,
int *cmatind,
double *cmatval,
double *lb,
double *ub);
int LP_change_bound(struct LP *lp,
int col,
char lower_or_upper,
double bnd);
int LP_optimize(struct LP *lp,
int *infeasible);
int LP_get_obj_val(struct LP *lp,
double *obj);
int LP_get_x(struct LP *lp,
double *x);
int LP_get_y(struct LP *lp,
double *y);
int LP_get_num_cols(const struct LP *lp);
int LP_get_num_rows(const struct LP *lp);
int LP_read_problem(struct LP *lp,
const char *filename);
int LP_relax(struct LP *lp);
int LP_get_column_types(struct LP *lp,
char *column_types);
int LP_get_tableau(struct LP *lp,
struct Row **rows,
int *cstat,
int *rstat,
double *ub,
double *lb);
int LP_unflip_row_coefficients(int *cstat,
double *lb,
double *ub,
struct Row *cut);
void LP_print_row(struct Row *row);
int LP_flip_row_coefficients(int *cstat,
double *lb,
double *ub,
struct Row *cut);
int LP_read_solution(struct LP *lp,
char *filename,
double *x);
void LP_disable_presolve(struct LP *lp);
int LP_write_solution(struct LP *lp,
char *filename);
int LP_write_basis(struct LP *lp,
char *filename);
int LP_read_basis(struct LP *lp,
char *filename);
int LP_write_sage_file(struct LP *lp,
double *x,
char *filename);
int LP_change_rhs(struct LP *lp, int index, double value);
#endif

26
multirow/include/multirow/mir.h
Unescape View File

@ -0,0 +1,26 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_MIR_H
#define MULTIROW_MIR_H
#include "lp.h"
int MIR_generate_cut(const struct Row *row,
const char *column_types,
struct Row *cut);
#endif //MULTIROW_MIR_H

90
multirow/include/multirow/params.h
Unescape View File

@ -0,0 +1,90 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PROJECT_PARAMS_H
#define PROJECT_PARAMS_H
/*
* Error margin for floating point comparisons.
*/
#define EPSILON 1e-5
/*
* Available log levels, in decreasing level of verboseness, are:
*
* LOG_LEVEL_VERBOSE
* LOG_LEVEL_DEBUG
* LOG_LEVEL_INFO
* LOG_LEVEL_WARNING
* LOG_LEVEL_ERROR
*/
#define LOG_LEVEL LOG_LEVEL_INFO
/*
* Maximum bounding-box size for naive algorithm
*/
#define MAX_BOX_SIZE 10000
/*
* Maximum number of sets that should be considered
*/
#define MAX_N_SETS 1000
/*
* Number of rays that should be generated per set.
*/
#define N_RAYS 100
#define ONLY_CUT -1
/*
* Time limit for the computation (user time, in seconds).
*/
#define MAX_TOTAL_TIME -1
/*
* Time limit for each CPLEX LP or MIP (in seconds).
*/
#define CPLEX_TIMEOUT 1.0
#define CG_TIMEOUT 900
#define MAX_N_RAYS 100
#define MAX_SELECTED_COMBINATIONS 5000
#define MAX_SELECTED_ROWS 300
#define MAX_LATTICE_POINTS 100000
/*
* Greedy cut parameters
*/
#define GREEDY_BIG_E 1e3
#define GREEDY_MAX_GAP 1e-4
#define MAX_CUT_DYNAMISM 1e8
#define INTEGRALITY_THRESHOLD 0.49
extern int BOOST_VAR;
extern double BOOST_FACTOR;
extern int DUMP_CUT;
extern int DUMP_CUT_N;
extern int ENABLE_LIFTING;
extern int MIN_N_ROWS;
extern int MAX_N_ROWS;
#define ERR_NO_CUT 2
#define ERR_MIP_TIMEOUT 3
#endif //PROJECT_PARAMS_H

34
multirow/include/multirow/rational.h
Unescape View File

@ -0,0 +1,34 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_RATIONAL_H
#define MULTIROW_RATIONAL_H
struct _rational
{
long num;
unsigned long den;
};
typedef struct _rational Rational[1];
long gcd(long a,
long b);
long lcm(long a,
long b);
#endif //MULTIROW_RATIONAL_H

33
multirow/include/multirow/stats.h
Unescape View File

@ -0,0 +1,33 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MULTIROW_STATS_H
#define MULTIROW_STATS_H
void STATS_init();
void STATS_set_input_filename(char *filename);
void STATS_set_obj_value(double obj);
void STATS_increment_generated_cuts();
void STATS_increment_added_cuts();
int STATS_print_yaml(char *filename);
void STATS_finish_round();
#endif

98
multirow/include/multirow/util.h
Unescape View File

@ -0,0 +1,98 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef _PROJECT_UTIL_H_
#define _PROJECT_UTIL_H_
#include <string.h>
#include <stdio.h>
#include "params.h"
void time_printf(const char *fmt,
...);
#define LOG_LEVEL_ERROR 10
#define LOG_LEVEL_WARNING 20
#define LOG_LEVEL_INFO 30
#define LOG_LEVEL_DEBUG 40
#define LOG_LEVEL_VERBOSE 50
#define if_error_level if(LOG_LEVEL >= LOG_LEVEL_ERROR)
#define if_warning_level if(LOG_LEVEL >= LOG_LEVEL_WARNING)
#define if_info_level if(LOG_LEVEL < LOG_LEVEL_INFO)
#define if_debug_level if(LOG_LEVEL >= LOG_LEVEL_DEBUG)
#define if_verbose_level if(LOG_LEVEL >= LOG_LEVEL_VERBOSE)
#if LOG_LEVEL < LOG_LEVEL_VERBOSE
#define log_verbose(...)
#else
#define log_verbose(...) time_printf( __VA_ARGS__)
#endif
#if LOG_LEVEL < LOG_LEVEL_DEBUG
#define log_debug(...)
#else
#define log_debug(...) time_printf( __VA_ARGS__)
#endif
#if LOG_LEVEL < LOG_LEVEL_INFO
#define log_info(...)
#else
#define log_info(...) time_printf( __VA_ARGS__)
#endif
#if LOG_LEVEL < LOG_LEVEL_WARNING
#define log_warn(...)
#else
#define log_warn(...) time_printf( __VA_ARGS__)
#endif
#if LOG_LEVEL < LOG_LEVEL_ERROR
#define log_error(...)
#else
#define log_error(...) time_printf( __VA_ARGS__)
#endif
#define abort_if(cond, msg) if(cond) { \
time_printf(msg " (%d) (%s:%d)\n", rval, __FILE__, __LINE__); \
rval = 1; goto CLEANUP; }
#define abort_iff(cond, msg, ...) if(cond) { \
time_printf(msg " (%d) (%s:%d)\n", __VA_ARGS__, rval, __FILE__, \
__LINE__); \
rval = 1; goto CLEANUP; }
#define UNUSED(x) (void)(x)
#define min(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)
double get_user_time(void);
double get_real_time();
double frac(double x);
extern FILE *LOG_FILE;
void progress_reset();
void progress_set_total(long total);
void progress_increment();
void progress_print();
void progress_title(char *title);
#endif

1021
multirow/src/cg.c
View File

File diff suppressed because it is too large Load Diff

134
multirow/src/double.c
Unescape View File

@ -0,0 +1,134 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <multirow/double.h>
#include <multirow/util.h>
int DOUBLE_sgn(double a)
{
return (a > 0 ? 1 : (a < 0 ? -1 : 0));
}
int DOUBLE_cmp(double a,
double b)
{
if(isnan(a) || isnan(b)) return 1;
double size = min(fabs(a) + fabs(b), DBL_MAX);
double zero_cutoff = EPSILON * size;
double diff = (a-b);
// Do not use relative epsilon for numbers close to zero
if(fabs(a) < EPSILON || fabs(b) < EPSILON) zero_cutoff = EPSILON;
return (fabs(diff) < zero_cutoff ? 0 : DOUBLE_sgn(diff));
}
int DOUBLE_geq(double a,
double b)
{
return DOUBLE_cmp(a, b) >= 0;
}
int DOUBLE_leq(double a,
double b)
{
return DOUBLE_cmp(a, b) <= 0;
}
int DOUBLE_eq(double a,
double b)
{
return DOUBLE_cmp(a, b) == 0;
}
int DOUBLE_neq(double a,
double b)
{
return DOUBLE_cmp(a, b) != 0;
}
double DOUBLE_max(double a,
double b)
{
return (DOUBLE_cmp(a,b) > 0 ? a : b);
}
double DOUBLE_random(double min, double max)
{
return min + ((double) rand() / ((double) RAND_MAX / (max - min)));
}
int DOUBLE_to_rational(double a,
long max_den,
Rational r)
{
int rval = 0;
abort_if(!isfinite(a), "a must be finite");
abort_if(max_den < 1, "max_den must be positive");
long shift = floor(a);
a = frac(a);
long n1 = 0, d1 = 1;
long n2 = 1, d2 = 1;
log_verbose("a=%.8lf\n", a);
while(d1 + d2 < max_den)
{
double mid = ((double) (n1 + n2) / (d1 + d2));
log_verbose(" d1=%ld %d\n", n1, d1);
log_verbose(" d2=%ld %d\n", n2, d2);
log_verbose(" mid=%.8lf\n", mid);
if(a > mid)
{
n1 = n1 + n2;
d1 = d1 + d2;
}
else
{
n2 = n1 + n2;
d2 = d1 + d2;
}
}
double err1 = fabs(a - ((double) n1 / d1));
double err2 = fabs(a - ((double) n2 / d2));
if(err1 < err2)
{
r->num = n1;
r->den = d1;
}
else
{
r->num = n2;
r->den = d2;
}
r->num += ((r->den) * shift);
CLEANUP:
return rval;
}

190
multirow/src/geometry.c
Unescape View File

@ -0,0 +1,190 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define _GNU_SOURCE
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <multirow/geometry.h>
#include <multirow/double.h>
#include <multirow/util.h>
struct SortPair
{
int index;
void *data;
};
#if defined(__APPLE__)
static int _qsort_cmp_angle_distance(void *p_orig,
const void *p1,
const void *p2)
#else
static int _qsort_cmp_angle_distance(const void *p1,
const void *p2,
void *p_orig)
#endif
{
double *orig = (double *) p_orig;
double *r1 = (double*) p1;
double *r2 = (double*) p2;
double norm1 = fabs(r1[0] - orig[0]) + fabs(r1[1] - orig[1]);
double norm2 = fabs(r2[0] - orig[0]) + fabs(r2[1] - orig[1]);
double angle1 = atan2(r1[0] - orig[0], r1[1] - orig[1]);
double angle2 = atan2(r2[0] - orig[0], r2[1] - orig[1]);
if(DOUBLE_neq(angle1, angle2))
return sign(angle2 - angle1);
else
return sign(norm1 - norm2);
}
static double turn_direction(const double *r1,
const double *r2,
const double *r3)
{
return (r2[0]-r1[0]) * (r3[1]-r1[1]) - (r2[1]-r1[1]) * (r3[0]-r1[0]);
}
int chull_2d(const double *original_points,
int npoints,
double *vertices,
int *nvertices)
{
int rval = 0;
int pivot_idx;
double pivot[2];
int *stack = 0;
int stack_length = 0;
double *points = 0;
stack = (int*) malloc(npoints * sizeof(int));
points = (double*) malloc(2 * npoints * sizeof(double));
abort_if(!stack, "could not allocate stack");
abort_if(!points, "could not allocate points");
memcpy(points, original_points, 2 * npoints * sizeof(double));
// find point with lowest y-coordinate
pivot_idx = 0;
memcpy(pivot, points, 2 * sizeof(double));
for (int i = 0; i < npoints; i++)
{
double *point = &points[2*i];
if(point[1] > pivot[1]) continue;
if(point[1] == pivot[1] && point[0] > pivot[0]) continue;
memcpy(pivot, point, 2 * sizeof(double));
pivot_idx = i;
}
// remove pivot
swap(points[2*(npoints-1)], points[2*pivot_idx], double);
swap(points[2*(npoints-1)+1], points[2*pivot_idx+1], double);
npoints--;
// sort points counterclockwise
#if defined(__APPLE__)
qsort_r(points, npoints, 2 * sizeof(double), pivot, _qsort_cmp_angle_distance);
#else
qsort_r(points, npoints, 2 * sizeof(double), _qsort_cmp_angle_distance,
pivot);
#endif
// adds pivot and first point to the stack
stack[stack_length++] = npoints;
stack[stack_length++] = 0;
int k = 1;
while(k < npoints)
{
double *p0 = &points[2*stack[stack_length-2]];
double *p1 = &points[2*stack[stack_length-1]];
double *p2 = &points[2*k];
double t = turn_direction(p0, p1, p2);
if(DOUBLE_leq(t, 0))
stack_length--;
if(DOUBLE_geq(t, 0))
{
stack[stack_length++] = k;
k++;
}
}
for(int i = 0; i < stack_length; i++)
{
int j = stack_length - i - 1;
vertices[2*j] = points[2*stack[i]];
vertices[2*j+1] = points[2*stack[i]+1];
}
*nvertices = stack_length;
CLEANUP:
if(points) free(points);
if(stack) free(stack);
return rval;
}
/*
* Find lambda such that lambda * p lies on the line connecting ray1 to ray2
*/
int scale_vector_to_line(const double *ray1,
const double *ray2,
const double *p,
double *lambda)
{
int rval = 0;
double a = ray1[0], b = ray1[1];
double c = ray2[0], d = ray2[1];
double norm1 = fabs(ray1[0]) + fabs(ray1[1]);
double norm2 = fabs(ray2[0]) + fabs(ray2[1]);
log_verbose("r1=%12.8e %12.8e\n", ray1[0], ray1[1]);
log_verbose("r2=%12.8e %12.8e\n", ray2[0], ray2[1]);
log_verbose("p=%12.8e %12.8e\n\n", p[0], p[1]);
if (norm1 > norm2)
{
swap(a, c, double);
swap(b, d, double);
swap(norm1, norm2, double);
}
if (DOUBLE_iszero(norm1 / norm2))
{
*lambda = (b * c - a * d) / (c * p[1] - d * p[0]);
}
else
{
*lambda = b * c - a * d;
*lambda /= p[0] * (b - d) - p[1] * (a - c);
}
CLEANUP:
return rval;
}

514
multirow/src/lfree2d.c
Unescape View File

@ -0,0 +1,514 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define _XOPEN_SOURCE 700
#include <math.h>
#include <stdio.h>
#include <multirow/double.h>
#include <multirow/util.h>
#include <multirow/lp.h>
#include <multirow/lfree2d.h>
int LFREE_2D_init(struct LFreeSet2D *set,
int n_vertices,
int n_lattice_points,
int n_halfspaces)
{
int rval = 0;
set->vertices = (double *) malloc(2 * n_vertices * sizeof(double));
set->lattice_points = (double *) malloc( 2 * n_lattice_points * sizeof(double));
set->halfspaces = (double *) malloc(2 * n_halfspaces * sizeof(double));
abort_if(!set->vertices, "could not allocate set->vertices");
abort_if(!set->lattice_points, "could not allocate set->lattice_points");
abort_if(!set->halfspaces, "could not allocate set->halfspaces");
CLEANUP:
return rval;
}
void LFREE_2D_free(struct LFreeSet2D *set)
{
if(set->vertices) free(set->vertices);
if(set->lattice_points) free(set->lattice_points);
if(set->halfspaces) free(set->halfspaces);
}
/**
* Computes the bounding box around the given set.
*
* @param[in] set the set whose bounding box should be computed
* @param[out] lb two-dimensional vector representing lower-left corner
* of the bounding box
* @param[out] ub two-dimensional vector representing upper-right corner
* of the bounding box
*/
int LFREE_2D_get_bounding_box(const struct LFreeSet2D *set,
int *lb,
int *ub)
{
int rval = 0;
ub[0] = ub[1] = INT_MIN;
lb[0] = lb[1] = INT_MAX;
for (int i = 0; i < set->n_vertices; i++)
{
double *vertex = &set->vertices[2 * i];
ub[0] = (int) fmax(ub[0], ceil(vertex[0]));
ub[1] = (int) fmax(ub[1], ceil(vertex[1]));
lb[0] = (int) fmin(lb[0], floor(vertex[0]));
lb[1] = (int) fmin(lb[1], floor(vertex[1]));
}
CLEANUP:
return rval;
}
int LFREE_2D_print_set(const struct LFreeSet2D *set)
{
int rval = 0;
log_debug(" f=%12.6lf %12.6lf\n", set->f[0], set->f[1]);
for (int i = 0; i < set->n_vertices; i++)
{
double *vertex = &set->vertices[2 * i];
log_debug(" v%-3d=%12.6lf %12.6lf\n", i, vertex[0], vertex[1]);
}
for (int i = 0; i < set->n_lattice_points; i++)
{
double *lattice_point = &set->lattice_points[2 * i];
log_debug(" t%-3d=%12.6lf %12.6lf\n", i, lattice_point[0], lattice_point[1]);
}
for(int i = 0; i < set->n_halfspaces; i++)
{
double *halfspace = &set->halfspaces[2 * i];
log_debug(" h%-3d=%12.6lf %12.6lf\n", i, halfspace[0], halfspace[1]);
}
CLEANUP:
return rval;
}
/**
* Computes the halfspace representation for a given set.
*
* To construct a LFreeSet2D, it is only necessary to specify
* the vertices and the interior point. This function computes
* the halfspace representation. It assumes that the vertices are
* sorted in either clockwise or counter-clockwise order.
* The set is modified in-place.
*
* @param set the set whose halfspace representation should
* be computed
*
*/
int LFREE_2D_compute_halfspaces(struct LFreeSet2D *set)
{
int rval = 0;
int k = 0;
int expected_sgn = 0;
double *f = set->f;
for (int i0 = 0; i0 < set->n_vertices; i0++)
{
int i1 = (i0 + 1) % set->n_vertices;
int i2 = (i0 + 2) % set->n_vertices;
double *v0 = &set->vertices[2 * i0];
double *v1 = &set->vertices[2 * i1];
double *v2 = &set->vertices[2 * i2];
double x0 = v0[0], y0 = v0[1];
double x1 = v1[0], y1 = v1[1];
double x2 = v2[0], y2 = v2[1];
int sgn = DOUBLE_sgn(
-x1 * y0 + x2 * y0 + x0 * y1 - x2 * y1 - x0 * y2 + x1 * y2);
abort_if(sgn == 0, "vertices should not be aligned");
if(expected_sgn == 0) expected_sgn = sgn;
abort_if(expected_sgn != sgn, "wrong orientation of vertices");
double *halfspace = &set->halfspaces[2 * (k++)];
halfspace[0] = y0 - y1;
halfspace[1] = x1 - x0;
double rhs = y0 * x1 - x0 * y1;
rhs -= (y0 - y1) * f[0];
rhs -= (x1 - x0) * f[1];
halfspace[0] /= rhs;
halfspace[1] /= rhs;
}
set->n_halfspaces = k;
CLEANUP:
return rval;
}
/**
* Translates a given set by a certain amount.
*
* @param set the set to be translated
* @param dx offset to be applied to the first component
* @param dy offset to be applied to the second component
*/
int LFREE_2D_translate_set(struct LFreeSet2D *set, double dx, double dy)
{
set->f[0] += dx;
set->f[1] += dy;
for (int i = 0; i < set->n_vertices; i++)
{
double *vertex = &set->vertices[2 * i];
vertex[0] += dx;
vertex[1] += dy;
}
for (int i = 0; i < set->n_lattice_points; i++)
{
double *lattice_point = &set->lattice_points[2 * i];
lattice_point[0] += dx;
lattice_point[1] += dy;
}
return 0;
}
/**
* Computes the inverse of the given matrix and left-multiplies it
* by the given ray. The matrix is given by [[x1,y1],[x2,y2]]. The
* ray is modified in-place.
*/
static void apply_inverse(double x1,
double y1,
double x2,
double y2,
double *ray)
{
double a = ray[0];
double b = ray[1];
double det = x2 * y1 - x1 * y2;
ray[0] = (b * x2 - a * y2) / det;
ray[1] = (a * y1 - b * x1) / det;
}
/**
* Applies a unimodular affine transformation to the set, making
* lattice points 1 and 2 orthogonal. Returns the transformation
* applied.
*
* @param[in] set the set to be transformed
* @param[out] m two-by-two matrix representing the transformation
* applied
*/
static int make_t1_t2_orthogonal(struct LFreeSet2D *set, double *m)
{
int rval = 0;
double x1 = set->lattice_points[2];
double y1 = set->lattice_points[3];
double x2 = set->lattice_points[4];
double y2 = set->lattice_points[5];
double det = x2 * y1 - x1 * y2;
abort_if(DOUBLE_iszero(det), "determinant should not be zero");
apply_inverse(x1, y1, x2, y2, set->f);
for (int i = 0; i < set->n_vertices; i++)
{
double *vertex = &set->vertices[2 * i];
apply_inverse(x1, y1, x2, y2, vertex);
}
for (int i = 0; i < set->n_lattice_points; i++)
{
double *lattice_point = &set->lattice_points[2 * i];
apply_inverse(x1, y1, x2, y2, lattice_point);
}
apply_inverse(x1, y1, x2, y2, &m[0]);
apply_inverse(x1, y1, x2, y2, &m[2]);
CLEANUP:
return rval;
}
static int apply_m_to_vect(double *v, const double *m)
{
double v0 = v[0], v1 = v[1];
v[0] = m[0] * v0 + m[1] * v1;
v[1] = m[2] * v0 + m[3] * v1;
return 0;
}
/**
* Applies a unimodular affine transformation to the set, making
* the four lattice points the unit square.
*
* Must be called on quadrilaterals only, and assumes that the three
* first lattice points already correspond to (0,0), (0,1) and (1,0).
* The set is modified in place. The transformation is also applied
* to a given matrix pre_m.
*/
static int make_square(struct LFreeSet2D *set, double *pre_m)
{
int rval = 0;
double m[4];
abort_if(set->n_vertices != 4, "make_square requires quadrilaterals");
double *t = &set->lattice_points[2 * 3];
if(t[0] == 1 && t[1] == 1)
{
m[0] = 1; m[1] = 0;
m[2] = 0; m[3] = 1;
}
else if(t[0] == 1 && t[1] == -1)
{
m[0] = 1; m[1] = 0;
m[2] = 1; m[3] = 1;
}
else if(t[0] == -1 && t[1] == 1)
{
m[0] = 1; m[1] = 1;
m[2] = 0; m[3] = 1;
}
else
{
abort_if(1, "invalid state");
}
rval = LFREE_2D_transform_set(set, m);
abort_if(rval, "LFREE_2D_transform_set failed");
apply_m_to_vect(&pre_m[0], m);
apply_m_to_vect(&pre_m[2], m);
CLEANUP:
return rval;
}
int find_thin_direction(const struct LFreeSet2D *set,
double *thin_direction,
double *width)
{
int rval = 0;
int n_candidates = 3;
double candidates[] = {
0, 1,
1, 0,
1, 1
};
double best_width = INFINITY;
double *best_candidate = 0;
for (int i = 0; i < n_candidates; i++)
{
double *d = &candidates[2 * i];
double max_val = -INFINITY;
double min_val = INFINITY;
for (int j = 0; j < set->n_vertices; j++)
{
double *vertex = &set->vertices[2 * j];
double val = vertex[0] * d[0] + vertex[1] * d[1];
max_val = fmax(max_val, val);
min_val = fmin(min_val, val);
}
double delta = fabs(max_val - min_val);
log_verbose("d=%.3lf %.3lf width=%.3lf\n", d[0], d[1], delta);
if(delta < best_width)
{
best_width = delta;
best_candidate = d;
}
}
memcpy(thin_direction, best_candidate, 2 * sizeof(double));
*width = best_width;
CLEANUP:
return rval;
}
int LFREE_2D_transform_set(struct LFreeSet2D *set, const double *m)
{
int rval = 0;
rval = apply_m_to_vect(set->f, m);
abort_if(rval, "apply_m_to_vect failed");
for (int i = 0; i < set->n_vertices; i++)
{
double *vertex = &set->vertices[2 * i];
rval = apply_m_to_vect(vertex, m);
abort_if(rval, "apply_m_to_vect failed");
}
for (int i = 0; i < set->n_lattice_points; i++)
{
double *lattice_point = &set->lattice_points[2 * i];
rval = apply_m_to_vect(lattice_point, m);
abort_if(rval, "apply_m_to_vect failed");
}
CLEANUP:
return rval;
}
int LFREE_2D_preprocess(struct LFreeSet2D *set, double *pre_m, double *center)
{
int rval = 0;
double direction[2], m[4], width;
double translate[2];
pre_m[0] = 1; pre_m[1] = 0;
pre_m[2] = 0; pre_m[3] = 1;
rval = LFREE_2D_translate_set(set, -set->lattice_points[0],
-set->lattice_points[1]);
abort_if(rval, "LFREE_2D_translate set failed");
rval = make_t1_t2_orthogonal(set, pre_m);
abort_if(rval, "make_t1_t2_orthogonal failed");
if(set->n_vertices == 4)
{
rval = make_square(set, pre_m);
abort_if(rval, "make_square failed");
}
rval = find_thin_direction(set, direction, &width);
abort_if(rval, "find_thin_direction failed");
log_debug("d=%.2lf %.2lf\n", direction[0], direction[1]);
if(direction[0] == 0 && direction[1] == 1)
{
m[0] = 1; m[1] = 0;
m[2] = 0; m[3] = 1;
translate[0] = translate[1] = 0;
}
else if(direction[0] == 1 && direction[1] == 0)
{
m[0] = 0; m[1] = 1;
m[2] = 1; m[3] = 0;
translate[0] = translate[1] = 0;
}
else if(direction[0] == 1 && direction[1] == 1)
{
m[0] = 1; m[1] = 0;
m[2] = -1; m[3] = -1;
translate[0] = 0;
translate[1] = 1;
}
else
{
abort_if(1, "invalid direction");
}
apply_m_to_vect(&pre_m[0], m);
apply_m_to_vect(&pre_m[2], m);
rval = LFREE_2D_transform_set(set, m);
abort_if(rval, "LFREE_2D_transform_set failed");
rval = LFREE_2D_translate_set(set, translate[0], translate[1]);
abort_if(rval, "LFREE_2D_translate_set failed");
double max_v1 = -INFINITY;
double min_v1 = INFINITY;
for(int i = 0; i < set->n_vertices; i++)
{
double *v = &set->vertices[2 * i];
max_v1 = fmax(max_v1, v[1]);
min_v1 = fmin(min_v1, v[1]);
}
*center = (max_v1 + min_v1) / 2;
width = max_v1 - min_v1;
log_debug("width = %.2lf\n", width);
abort_iff(width >= 3, "lattice width too large: %.2lf", width);
CLEANUP:
return rval;
}
int LFREE_2D_read_next(FILE *file, struct LFreeSet2D *set)
{
int rval = 0;
ssize_t count;
double *f = set->f;
count = fscanf(file, "%lf %lf ", &f[0], &f[1]);
abort_if(count != 2, "could not read f");
count = fscanf(file, "%d ", &set->n_vertices);
abort_if(count != 1, "could not read n_vertices");
abort_if(set->n_vertices != 3 && set->n_vertices != 4,
"only triangles or quadrilaterals supported");
for(int i = 0; i < set->n_vertices; i++)
{
double *vertex = &set->vertices[2 * i];
count = fscanf(file, "%lf %lf ", &vertex[0], &vertex[1]);
abort_iff(count != 2, "could not read vertex %d", i+1);
}
count = fscanf(file, "%d ", &set->n_lattice_points);
abort_if(count != 1, "could not read n_lattice_points");
abort_if(set->n_lattice_points != set->n_vertices, "number of vertices "
"should be the same as number of lattice points");
for(int i = 0; i < set->n_lattice_points; i++)
{
double *lattice_point = &set->lattice_points[2 * i];
count = fscanf(file, "%lf %lf ", &lattice_point[0], &lattice_point[1]);
abort_iff(count != 2, "could not read lattice_point %d", i+1);
}
CLEANUP:
return rval;
}

724
multirow/src/lp.c
Unescape View File

@ -0,0 +1,724 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ilcplex/cplex.h>
#include <multirow/lp.h>
#include <multirow/util.h>
#include <multirow/double.h>
int LP_open(struct LP *lp)
{
int rval = 0;
lp->cplex_lp = 0;
lp->cplex_env = CPXopenCPLEX(&rval);
abort_if(rval, "CPXopenCPLEX failed");
CPXsetintparam(lp->cplex_env, CPX_PARAM_DATACHECK, CPX_ON);
CPXsetintparam(lp->cplex_env, CPX_PARAM_NUMERICALEMPHASIS, CPX_ON);
CPXsetlogfile(lp->cplex_env, 0);
CLEANUP:
return rval;
}
void LP_disable_presolve(struct LP *lp)
{
CPXsetintparam(lp->cplex_env, CPX_PARAM_PREIND, CPX_OFF);
CPXsetintparam(lp->cplex_env, CPX_PARAM_REPEATPRESOLVE, CPX_OFF);
}
int LP_create(struct LP *lp,
const char *name)
{
int rval = 0;
char nambuf[MAX_NAME_LENGTH];
abort_if(!lp->cplex_env, "cplex_env is null");
strncpy(nambuf, name, MAX_NAME_LENGTH);
nambuf[MAX_NAME_LENGTH - 1] = '\0';
lp->cplex_lp = CPXcreateprob(lp->cplex_env, &rval, nambuf);
abort_if(rval, "CPXcreateprob failed");
CLEANUP:
return rval;
}
void LP_free(struct LP *lp)
{
if (!lp) return;
if (!lp->cplex_env) return;
if (lp->cplex_lp)
CPXfreeprob(lp->cplex_env, &(lp->cplex_lp));
CPXcloseCPLEX(&lp->cplex_env);
lp->cplex_env = 0;
}
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,
int newrows,
int newnz,
double *rhs,
char *sense,
int *rmatbeg,
int *rmatind,
double *rmatval)
{
int rval = 0;
rval = CPXaddrows(lp->cplex_env, lp->cplex_lp, 0, newrows, newnz, rhs,
sense, rmatbeg, rmatind, rmatval, 0, 0);
abort_if(rval, "CPXaddrows failed");
CLEANUP:
return rval;
}
int LP_add_row(struct LP *lp,
struct Row *row)
{
int rval = 0;
char sense = 'L';
int rmatbeg = 0;
rval = CPXaddrows(lp->cplex_env, lp->cplex_lp, 0, 1, row->nz, &row->pi_zero,
&sense, &rmatbeg, row->indices, row->pi, 0, 0);
abort_if(rval, "CPXaddrows failed");
CLEANUP:
return rval;
}
int LP_add_cols(struct LP *lp,
int newcols,
int newnz,
double *obj,
int *cmatbeg,
int *cmatind,
double *cmatval,
double *lb,
double *ub)
{
int rval = 0;
rval = CPXaddcols(lp->cplex_env, lp->cplex_lp, newcols, newnz, obj, cmatbeg,
cmatind, cmatval, lb, ub, 0);
abort_if(rval, "CPXaddcols failed");
CLEANUP:
return rval;
}
int LP_change_bound(struct LP *lp,
int col,
char lower_or_upper,
double bnd)
{
int rval = 0;
rval = CPXchgbds(lp->cplex_env, lp->cplex_lp, 1, &col, &lower_or_upper,
&bnd);
abort_if(rval, "CPXchgbds failed");
CLEANUP:
return rval;
}
int LP_get_obj_val(struct LP *lp,
double *obj)
{
int rval = 0;
rval = CPXgetobjval(lp->cplex_env, lp->cplex_lp, obj);
abort_if(rval, "CPXgetobjval failed");
CLEANUP:
return rval;
}
int LP_get_x(struct LP *lp,
double *x)
{
int rval = 0;
int ncols = CPXgetnumcols(lp->cplex_env, lp->cplex_lp);
abort_if(!ncols, "No columns in LP");
rval = CPXgetx(lp->cplex_env, lp->cplex_lp, x, 0, ncols - 1);
abort_if(rval, "CPXgetx failed");
CLEANUP:
return rval;
}
int LP_get_y(struct LP *lp,
double *y)
{
int rval = 0;
int nrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
abort_if(!nrows, "No rows in LP");
rval = CPXgetpi(lp->cplex_env, lp->cplex_lp, y, 0, nrows - 1);
abort_iff(rval, "CPXgetpi failed (errno = %d)", rval);
CLEANUP:
return rval;
}
int LP_get_num_cols(const struct LP *lp)
{
return CPXgetnumcols(lp->cplex_env, lp->cplex_lp);
}
int LP_get_num_rows(const struct LP *lp)
{
return CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
}
int LP_new_row(struct LP *lp,
char sense,
double rhs)
{
int rval = 0;
rval = CPXnewrows(lp->cplex_env, lp->cplex_lp, 1, &rhs, &sense, 0, 0);
abort_if(rval, "CPXnewrows failed");
CLEANUP:
return rval;
}
int LP_new_col(struct LP *lp,
double obj,
double lb,
double ub,
char type)
{
int rval = 0;
rval = CPXnewcols(lp->cplex_env, lp->cplex_lp, 1, &obj, &lb, &ub, &type, 0);
abort_if(rval, "CPXnewcols failed");
CLEANUP:
return rval;
}
int LP_optimize(struct LP *lp,
int *infeasible)
{
int rval = 0, solstat;
*infeasible = 0;
// LP_SOLVE_COUNT++;
if_verbose_level
{
int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
int numcols = CPXgetnumcols(lp->cplex_env, lp->cplex_lp);
time_printf("Optimizing MILP (%d rows %d cols)...\n", numrows, numcols);
}
double initial_time = get_user_time();
int problem_type = CPXgetprobtype(lp->cplex_env, lp->cplex_lp);
switch (problem_type)
{
case CPXPROB_LP:
rval = CPXdualopt(lp->cplex_env, lp->cplex_lp);
log_verbose(" dual opt\n");
abort_if(rval, "CPXdualopt failed");
break;
case CPXPROB_MILP:
case CPXPROB_FIXEDMILP:
rval = CPXmipopt(lp->cplex_env, lp->cplex_lp);
log_verbose(" mip opt\n");
abort_if(rval, "CPXmipopt failed");
break;
default:
abort_iff(1, "Invalid problem type: %d", problem_type);
}
solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp);
switch (solstat)
{
case CPX_STAT_INFEASIBLE:
case CPXMIP_INFEASIBLE:
log_verbose(" infeasible\n");
*infeasible = 1;
goto CLEANUP;
case CPXMIP_OPTIMAL:
case CPXMIP_OPTIMAL_INFEAS:
case CPXMIP_OPTIMAL_TOL:
case CPX_STAT_OPTIMAL:
case CPX_STAT_OPTIMAL_INFEAS:
case CPX_STAT_UNBOUNDED:
break;
default:
abort_iff(1, "Invalid solution status: %d", solstat);
}
double objval;
rval = LP_get_obj_val(lp, &objval);
abort_if(rval, "LP_get_obj_val failed");
log_verbose(" obj val = %.4lf\n", objval);
log_verbose(" time = %.4lf\n", get_user_time() - initial_time);
CLEANUP:
return rval;
}
int LP_write(struct LP *lp,
const char *fname)
{
int rval = 0;
char nambuf[MAX_NAME_LENGTH];
FILE *f = fopen(fname, "w");
abort_iff(!f, "could not open file %s", fname);
fclose(f);
strncpy(nambuf, fname, MAX_NAME_LENGTH);
nambuf[MAX_NAME_LENGTH - 1] = '\0';
log_info("Writing LP to file %s...\n", fname);
rval = CPXwriteprob(lp->cplex_env, lp->cplex_lp, nambuf, "RLP");
abort_if(rval, "CPXwriteprob failed");
CLEANUP:
return rval;
}
int LP_read_problem(struct LP *lp,
const char *filename)
{
int rval = 0;
log_info("Reading problem %s...\n", filename);
rval = CPXreadcopyprob(lp->cplex_env, lp->cplex_lp, filename, 0);
abort_if(rval, "CPXreadcopyprob failed");
CLEANUP:
return rval;
}
int LP_relax(struct LP *lp)
{
int rval = 0;
rval = CPXchgprobtype(lp->cplex_env, lp->cplex_lp, CPXPROB_LP);
abort_if(rval, "CPXchgprobtype failed");
CLEANUP:
return rval;
}
int LP_get_column_types(struct LP *lp,
char *column_types)
{
int rval = 0;
char *cplex_ctype = 0;
int ncols = LP_get_num_cols(lp);
cplex_ctype = (char *) malloc(ncols * sizeof(char));
abort_if(!cplex_ctype, "could not allocate cplex_ctype");
rval = CPXgetctype(lp->cplex_env, lp->cplex_lp, cplex_ctype, 0, ncols - 1);
abort_if(rval, "CPXgetctype failed");
for (int i = 0; i < ncols; i++)
{
switch (cplex_ctype[i])
{
case CPX_BINARY:
case CPX_INTEGER:
column_types[i] = MILP_INTEGER;
break;
case CPX_CONTINUOUS:
column_types[i] = MILP_CONTINUOUS;
break;
default:
abort_iff(1, "Invalid column type: %d", cplex_ctype[i]);
}
}
CLEANUP:
if (cplex_ctype) free(cplex_ctype);
return rval;
}
int LP_get_tableau(struct LP *lp,
struct Row **rows,
int *cstat,
int *rstat,
double *ub,
double *lb)
{
int rval = 0;
int nrows = LP_get_num_rows(lp);
int ncols = LP_get_num_cols(lp);
int *head = 0;
double *rhs = 0;
double *pi = 0;
double *slacks = 0;
rhs = (double *) malloc(nrows * sizeof(double));
head = (int *) malloc(nrows * sizeof(int));
slacks = (double*) malloc(nrows * sizeof(double));
abort_if(!head, "could not allocate head");
abort_if(!rhs, "could not allocate rhs");
abort_if(!slacks, "could not allocate slacks");
rval = CPXgetbhead(lp->cplex_env, lp->cplex_lp, head, rhs);
abort_if(rval, "CPXgetbhead failed");
rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, nrows-1);
abort_if(rval, "CPXgetslack failed");
rval = CPXgetub(lp->cplex_env, lp->cplex_lp, ub, 0, ncols - 1);
abort_if(rval, "CPXgetub failed");
rval = CPXgetlb(lp->cplex_env, lp->cplex_lp, lb, 0, ncols - 1);
abort_if(rval, "CPXgetlb failed");
rval = CPXgetbase(lp->cplex_env, lp->cplex_lp, cstat, rstat);
abort_if(rval, "CPXgetbase failed");
for (int i = 0; i < nrows; i++)
abort_if(!DOUBLE_iszero(slacks[i]), "model contains slack variables");
pi = (double *) malloc(ncols * sizeof(double));
abort_if(!pi, "could not allocate pi");
for (int i = 0; i < nrows; i++)
{
int nz = 0;
rows[i] = (struct Row *) malloc(sizeof(struct Row));
abort_if(!rows[i], "could not allocate rows[i]");
struct Row *row = rows[i];
rval = CPXbinvarow(lp->cplex_env, lp->cplex_lp, i, pi);
abort_if(rval, "CPXbinvarow failed");
for (int j = 0; j < ncols; j++)
{
if (fabs(pi[j]) < EPSILON)
continue;
if (cstat[j] == CPX_AT_LOWER)
rhs[i] += lb[j] * pi[j];
if (cstat[j] == CPX_AT_UPPER)
rhs[i] += ub[j] * pi[j];
nz++;
}
row->nz = nz;
row->pi_zero = rhs[i];
row->head = head[i];
row->pi = (double *) malloc(nz * sizeof(double));
row->indices = (int *) malloc(nz * sizeof(int));
abort_if(!row->pi, "could not allocate row->pi");
abort_if(!row->indices, "could not allocate row->indices");
if (fabs(row->pi_zero) < EPSILON)
row->pi_zero = 0;
int k = 0;
for (int j = 0; j < ncols; j++)
{
if (fabs(pi[j]) < EPSILON)
continue;
row->pi[k] = pi[j];
row->indices[k++] = j;
}
rval = LP_flip_row_coefficients(cstat, lb, ub, row);
abort_if(rval, "LP_flip_row_coefficients failed");
}
CLEANUP:
if(slacks) free(slacks);
if (head) free(head);
if (rhs) free(rhs);
if (pi) free(pi);
return rval;
}
int LP_flip_row_coefficients(int *cstat,
double *lb,
double *ub,
struct Row *cut)
{
for (int j = 0; j < cut->nz; j++)
{
int idx = cut->indices[j];
double pij = cut->pi[j];
if (cstat[idx] == CPX_AT_LOWER)
cut->pi_zero -= lb[idx] * pij;
if (cstat[idx] == CPX_AT_UPPER)
{
cut->pi_zero -= ub[idx] * pij;
pij = -pij;
}
cut->indices[j] = idx;
cut->pi[j] = pij;
}
return 0;
}
int LP_unflip_row_coefficients(int *cstat,
double *lb,
double *ub,
struct Row *cut)
{
for (int j = 0; j < cut->nz; j++)
{
int idx = cut->indices[j];
double pij = cut->pi[j];
if (cstat[idx] == CPX_AT_LOWER)
cut->pi_zero += lb[idx] * pij;
if (cstat[idx] == CPX_AT_UPPER)
{
pij = -pij;
cut->pi_zero += ub[idx] * pij;
}
cut->indices[j] = idx;
cut->pi[j] = pij;
}
return 0;
}
void LP_print_row(struct Row *row)
{
for (int i = 0; i < row->nz; i++)
printf("%.2lfx%d ", row->pi[i], row->indices[i]);
printf(" <= %.2lf\n", row->pi_zero);
}
int LP_read_solution(struct LP *lp,
char *filename,
double *x)
{
int rval = 0;
int ncols = LP_get_num_cols(lp);
log_info("Reading solution from %s\n", filename);
FILE *fsol = fopen(filename, "rb");
abort_iff(!fsol, "Could not open file %s", filename)
for (int i = 0; i < ncols; i++)
{
int count = fscanf(fsol, "%le", &x[i]);
abort_iff(count != 1, "Unexpected EOF when reading %s", filename);
}
fclose(fsol);
CLEANUP:
return rval;
}
int LP_write_solution(struct LP *lp,
char *filename)
{
int rval = 0;
long ncols = LP_get_num_cols(lp);
double *x = 0;
x = (double *) malloc(ncols * sizeof(double));
abort_if(!x, "could not allocate x");
rval = LP_get_x(lp, x);
abort_if(rval, "LP_get_x failed");
log_info("Writing solution to file %s\n", filename);
FILE *fsol = fopen(filename, "w");
abort_iff(!fsol, "Could not open file %s", filename)
for (int i = 0; i < ncols; i++)
fprintf(fsol, "%.12e\n", x[i]);
fclose(fsol);
CLEANUP:
if (x) free(x);
return rval;
}
int LP_write_basis(struct LP *lp,
char *filename)
{
int rval = 0;
log_info("Writing basis to file %s...\n", filename);
rval = CPXmbasewrite(lp->cplex_env, lp->cplex_lp, filename);
abort_if(rval, "CPXmbasewrite failed");
CLEANUP:
return rval;
}
int LP_read_basis(struct LP *lp,
char *filename)
{
int rval = 0;
log_info("Reading basis from file %s...\n", filename);
rval = CPXreadcopybase(lp->cplex_env, lp->cplex_lp, filename);
abort_if(rval, "CPXreadcopybase failed");
CLEANUP:
return rval;
}
int LP_write_sage_file(struct LP *lp,
double *x,
char *filename)
{
int rval = 0;
FILE *fsage = 0;
int ncols = LP_get_num_cols(lp);
int nrows = LP_get_num_rows(lp);
int nz;
int dummy;
int *indices = 0;
double *rhs = 0;
double *values = 0;
int *head = 0;
fsage = fopen(filename, "w");
abort_if(!fsage, "fopen failed");
indices = (int*) malloc(ncols * sizeof(int));
values = (double*) malloc(ncols * sizeof(double));
rhs = (double*) malloc(nrows * sizeof(double));
head = (int*) malloc(nrows * sizeof(int));
abort_if(!indices, "could not allocate indices");
abort_if(!values, "could not allocate values");
abort_if(!rhs, "could not allocate rhs");
abort_if(!head, "could not allocate head");
rval = CPXgetrhs(lp->cplex_env, lp->cplex_lp, rhs, 0, nrows-1);
abort_if(rval, "CPXgetrhs failed");
fprintf(fsage, "b = vector([");
for (int i = 0; i < nrows; i++)
fprintf(fsage, " %.8lf,\n", rhs[i]);
fprintf(fsage, "])\n");
fprintf(fsage, "A = matrix(QQ, %d, %d, sparse=True)\n", nrows, ncols);
for (int i = 0; i < nrows; i++)
{
rval = CPXgetrows(lp->cplex_env, lp->cplex_lp, &nz, &dummy, indices,
values, ncols, &dummy, i, i);
abort_if(rval, "CPXgetrows failed");
for (int k = 0; k < nz; k++)
fprintf(fsage, "A[%5d,%5d] = %20.8lf\n", i, indices[k], values[k]);
}
rval = CPXgetbhead(lp->cplex_env, lp->cplex_lp, head, rhs);
abort_if(rval, "CPXgetbhead failed");
fprintf(fsage, "B=[");
for (int i = 0; i < nrows; i++)
fprintf(fsage, " %d,\n", head[i]);
fprintf(fsage, "]\n");
fprintf(fsage, "x=vector([");
for (int i = 0; i < ncols; i++)
fprintf(fsage, " %.8lf,\n", x[i]);
fprintf(fsage, "])\n");
fclose(fsage);
CLEANUP:
if(indices) free(indices);
if(values) free(values);
if(rhs) free(rhs);
if(head) free(head);
return rval;
}
int LP_change_rhs(struct LP *lp, int index, double value)
{
int rval = 0;
rval = CPXchgrhs(lp->cplex_env, lp->cplex_lp, 1, &index, &value);
abort_if(rval, "CPXchgrhs failed");
CLEANUP:
return rval;
}

63
multirow/src/mir.c
Unescape View File

@ -0,0 +1,63 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <multirow/lp.h>
#include <multirow/mir.h>
#include <multirow/util.h>
static double h(double a)
{
return fmax(a, 0);
}
static double f(double a,
double b)
{
return frac(b) * floor(a) + fmin(frac(a), frac(b));
}
int MIR_generate_cut(const struct Row *row,
const char *column_types,
struct Row *cut)
{
int rval = 0;
cut->pi = (double*) malloc(row->nz * sizeof(double));
cut->indices = (int*) malloc(row->nz * sizeof(int));
abort_if(!cut->pi, "could not allocate cut->pi");
abort_if(!cut->indices, "could not allocate cut->indices");
cut->head = 0;
cut->nz = row->nz;
cut->pi_zero = -frac(row->pi_zero) * ceil(row->pi_zero);
for (int i = 0; i < cut->nz; i++)
{
int idx = row->indices[i];
double value = row->pi[i];
cut->indices[i] = idx;
if (column_types[idx] == MILP_INTEGER)
cut->pi[i] = -f(value, row->pi_zero);
else
cut->pi[i] = -h(value);
}
CLEANUP:
return rval;
}

37
multirow/src/rational.c
Unescape View File

@ -0,0 +1,37 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
long gcd(long a,
long b)
{
long temp;
while (b != 0)
{
temp = a % b;
a = b;
b = temp;
}
return a;
}
long lcm(long a,
long b)
{
long g = gcd(a, b);
return (a/g) * b;
}

120
multirow/src/stats.c
Unescape View File

@ -0,0 +1,120 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define MAX_LENGTH 1000
#define MAX_ROUNDS 100
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <multirow/util.h>
char input_filename[MAX_LENGTH];
unsigned long long generated_cuts_count[MAX_LENGTH];
unsigned long long added_cuts_count[MAX_LENGTH];
int n_rounds = 0;
unsigned long long total_added_cuts = 0;
unsigned long long total_generated_cuts = 0;
double obj_value[MAX_ROUNDS];
double prev_time;
double runtime[MAX_ROUNDS];
void STATS_init()
{
for(int i = 0; i < MAX_ROUNDS; i++)
{
obj_value[i] = NAN;
generated_cuts_count[i] = 0;
added_cuts_count[i] = 0;
runtime[i] = 0;
}
prev_time = get_user_time();
}
void STATS_set_input_filename(char *filename)
{
strcpy(input_filename, filename);
}
void STATS_set_obj_value(double obj)
{
obj_value[n_rounds] = obj;
}
void STATS_increment_generated_cuts()
{
generated_cuts_count[n_rounds]++;
total_generated_cuts++;
}
void STATS_increment_added_cuts()
{
added_cuts_count[n_rounds]++;
total_added_cuts++;
}
void STATS_finish_round()
{
double now = get_user_time();
runtime[n_rounds] = now - prev_time;
prev_time = now;
n_rounds++;
}
int STATS_print_yaml(char *filename)
{
int rval = 0;
FILE *f = fopen(filename, "w");
abort_iff(!f, "could not open file %s", filename);
fprintf(f, "input-filename:\n %s\n", input_filename);
fprintf(f, "obj-value:\n");
for(int i = 0; i < n_rounds; i++)
fprintf(f, " %d: %.6lf\n", i, obj_value[i]);
fprintf(f, "generated-cuts:\n");
for(int i = 0; i < n_rounds; i++)
fprintf(f, " %d: %lld\n", i, generated_cuts_count[i]);
fprintf(f, "added-cuts:\n");
for(int i = 0; i < n_rounds; i++)
fprintf(f, " %d: %lld\n", i, added_cuts_count[i]);
fprintf(f, "user-cpu-time:\n");
for(int i = 0; i < n_rounds; i++)
fprintf(f, " %d: %.3lf\n", i, runtime[i]);
fprintf(f, "time-per-cut:\n");
for(int i = 0; i < n_rounds; i++)
{
if(generated_cuts_count[i] > 0)
fprintf(f, " %d: %.3lf\n", i, runtime[i] / generated_cuts_count[i]);
}
fclose(f);
CLEANUP:
return rval;
}

134
multirow/src/util.c
Unescape View File

@ -0,0 +1,134 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdio.h>
#include <sys/resource.h>
#include <sys/ioctl.h>
#include <stdarg.h>
#include <time.h>
#include <math.h>
#include <unistd.h>
#include <multirow/util.h>
FILE *LOG_FILE = 0;
double INITIAL_TIME = 0;
double get_user_time()
{
struct rusage ru;
getrusage(RUSAGE_SELF, &ru);
return ((double) ru.ru_utime.tv_sec)
+ ((double) ru.ru_utime.tv_usec) / 1000000.0;
}
double get_real_time()
{
return (double) time (0);
}
void time_printf(const char *fmt,
...)
{
if (INITIAL_TIME == 0)
INITIAL_TIME = get_user_time();
printf("[%10.2lf] ", get_user_time() - INITIAL_TIME);
va_list args;
va_start(args, fmt);
vprintf(fmt, args);
va_end(args);
fflush(stdout);
if(LOG_FILE)
{
fprintf(LOG_FILE, "[%10.2lf] ", get_user_time() - INITIAL_TIME);
va_list args;
va_start(args, fmt);
vfprintf(LOG_FILE, fmt, args);
va_end(args);
fflush(LOG_FILE);
}
}
double frac(double x)
{
return x - floor(x);
}
static long eta_count;
static long eta_total;
static time_t eta_start;
static time_t eta_last;
static char eta_title[1000] = {0};
void progress_reset()
{
eta_count = 0;
time(&eta_start);
eta_last = eta_start;
}
void progress_set_total(long total)
{
eta_total = total;
}
void progress_increment()
{
eta_count++;
}
void progress_title(char *title)
{
strncpy(eta_title, title, 1000);
}
void progress_print()
{
struct winsize w;
ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
if (eta_count == 0) return;
time_t eta_now;
time(&eta_now);
if(fabs(difftime(eta_now, eta_last)) < 30.0)
return;
eta_last = eta_now;
double diff = difftime(eta_now, eta_start);
double eta = diff / eta_count * eta_total;
int length = 10;
time_t eta_date = eta_start + eta;
struct tm *ttm = localtime(&eta_date);
fprintf(stdout,
"[%-40s] %*s%3.0f%% ETA: %04d-%02d-%02d %02d:%02d %*ld / %*ld\n",
eta_title, 0, "", 100.0 * eta_count / eta_total, ttm->tm_year +
1900, ttm->tm_mon + 1, ttm->tm_mday, ttm->tm_hour, ttm->tm_min,
length, eta_count, length, eta_total);
fflush(stdout);
}

267
multirow/tests/cg-test.cpp
Unescape View File

@ -0,0 +1,267 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#define TEST_SOURCE
#include <gtest/gtest.h>
#include <math.h>
extern "C" {
#include <multirow/cg.h>
#include "../src/cg.c"
}
int BOOST_VAR = -1;
double BOOST_FACTOR = 1.0;
TEST(CGTest, next_combination_test_1)
{
int rval = 0;
int v[] = { 1, 0 };
int n = 10;
int k = 2;
int finished;
int count = 0;
do
{
count++;
next_combination(n, k, v, 0, &finished);
}
while(!finished);
EXPECT_EQ(count, 45);
CLEANUP:
if(rval) FAIL();
}
TEST(CGTest, next_combination_test_2)
{
int rval = 0;
int v[] = { 2, 1, 0 };
int n = 5;
int k = 3;
int finished;
int count = 0;
do
{
count++;
next_combination(n, k, v, 0, &finished);
}
while(!finished);
EXPECT_EQ(count, 10);
CLEANUP:
if(rval) FAIL();
}
TEST(CGTest, check_rays_parallel_test)
{
int rval = 0;
double r1[] = { 1.0, 1.0, 1.0 };
double r2[] = { 2.0, 2.0, 2.0 };
double r3[] = { 0.0, 0.0, 0.0 };
double r4[] = { -1.0, -1.0, -1.0 };
int match;
double scale;
rval = check_rays_parallel(3, r1, r2, &match, &scale);
abort_if(rval, "check_rays_parallel failed");
EXPECT_TRUE(match);
EXPECT_DOUBLE_EQ(scale, 0.5);
rval = check_rays_parallel(3, r1, r3, &match, &scale);
abort_if(rval, "check_rays_parallel failed");
EXPECT_FALSE(match);
rval = check_rays_parallel(3, r1, r4, &match, &scale);
abort_if(rval, "check_rays_parallel failed");
EXPECT_FALSE(match);
CLEANUP:
if(rval) FAIL();
}
TEST(CGTest, extract_rays_from_rows_test)
{
int rval = 0;
double pi1[] = { 1.0, 1.0, 1.0, 2.0, 1.0 };
int indices1[] = { 1, 7, 8, 12, 14 };
struct Row row1 =
{
.nz = 5,
.head = 1,
.pi_zero = 1.0,
.pi = pi1,
.indices = indices1
};
double pi2[] = { 1.0, 1.0, 2.0, 1.0, 1.0, 1.0 };
int indices2[] = { 5, 8, 12, 13, 14, 16 };
struct Row row2 =
{
.nz = 6,
.head = 13,
.pi_zero = 1.0,
.pi = pi2,
.indices = indices2
};
double pi3[] = { 1.0, 1.0, 1.0, 1.0 };
int indices3[] = { 3, 7, 10, 16 };
struct Row row3 =
{
.nz = 4,
.head = 7,
.pi_zero = 1.0,
.pi = pi3,
.indices = indices3
};
struct Row *rows[] = { &row1, &row2, &row3 };
int nz;
int nrays;
int indices[1000];
int variable_to_ray[1000];
double rays[1000];
double ray_scale[1000];
rval = CG_extract_rays_from_rows(3, rows, rays, &nrays, variable_to_ray,
ray_scale, indices, &nz);
abort_if(rval, "CG_extract_rays_from_rows failed");
EXPECT_EQ(nrays, 4);
EXPECT_EQ(nz, 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[3], 0.0);
EXPECT_DOUBLE_EQ(rays[4], -1.0);
EXPECT_DOUBLE_EQ(rays[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[ 9], 0.0);
EXPECT_DOUBLE_EQ(rays[10], -1.0);
EXPECT_DOUBLE_EQ(rays[11], -1.0);
EXPECT_EQ(indices[0], 3);
EXPECT_EQ(indices[1], 5);
EXPECT_EQ(indices[2], 8);
EXPECT_EQ(indices[3], 10);
EXPECT_EQ(indices[4], 12);
EXPECT_EQ(indices[5], 14);
EXPECT_EQ(indices[6], 16);
EXPECT_EQ(ray_scale[0], 1.0);
EXPECT_EQ(ray_scale[1], 1.0);
EXPECT_EQ(ray_scale[2], 0.5);
EXPECT_EQ(ray_scale[3], 1.0);
EXPECT_EQ(ray_scale[4], 1.0);
EXPECT_EQ(ray_scale[5], 0.5);
EXPECT_EQ(ray_scale[6], 1.0);
EXPECT_EQ(variable_to_ray[0], 0);
EXPECT_EQ(variable_to_ray[1], 1);
EXPECT_EQ(variable_to_ray[2], 2);
EXPECT_EQ(variable_to_ray[3], 0);
EXPECT_EQ(variable_to_ray[4], 2);
EXPECT_EQ(variable_to_ray[5], 2);
EXPECT_EQ(variable_to_ray[6], 3);
CLEANUP:
if(rval) FAIL();
}
TEST(CGTest, boost_variable_test)
{
int rval = 0;
int nrows = 3;
double rays[] = {
0, 0, -1,
0, -1, 0,
-2, -2, 0,
0, -1, -1
};
int nz = 7;
int indices[] = { 3, 5, 8, 10, 12, 14, 16 };
int variable_to_ray[] = { 0, 1, 2, 0, 2, 2, 3 };
double ray_scale[] = { 1, 1, 0.5, 1, 1, 0.5, 1 };
rval = CG_boost_variable(8, 200, nrows, rays, variable_to_ray, ray_scale, indices, nz);
abort_if(rval, "boost_variable failed");
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[3], 0.0);
EXPECT_DOUBLE_EQ(rays[4], -1.0);
EXPECT_DOUBLE_EQ(rays[5], 0.0);
EXPECT_DOUBLE_EQ(rays[6], -400.0);
EXPECT_DOUBLE_EQ(rays[7], -400.0);
EXPECT_DOUBLE_EQ(rays[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_EQ(indices[0], 3);
EXPECT_EQ(indices[1], 5);
EXPECT_EQ(indices[2], 8);
EXPECT_EQ(indices[3], 10);
EXPECT_EQ(indices[4], 12);
EXPECT_EQ(indices[5], 14);
EXPECT_EQ(indices[6], 16);
EXPECT_EQ(ray_scale[0], 1.0);
EXPECT_EQ(ray_scale[1], 1.0);
EXPECT_EQ(ray_scale[2], 0.0025);
EXPECT_EQ(ray_scale[3], 1.0);
EXPECT_EQ(ray_scale[4], 0.005);
EXPECT_EQ(ray_scale[5], 0.0025);
EXPECT_EQ(ray_scale[6], 1.0);
EXPECT_EQ(variable_to_ray[0], 0);
EXPECT_EQ(variable_to_ray[1], 1);
EXPECT_EQ(variable_to_ray[2], 2);
EXPECT_EQ(variable_to_ray[3], 0);
EXPECT_EQ(variable_to_ray[4], 2);
EXPECT_EQ(variable_to_ray[5], 2);
EXPECT_EQ(variable_to_ray[6], 3);
CLEANUP:
if(rval) FAIL();
}

164
multirow/tests/double-test.cpp
Unescape View File

@ -0,0 +1,164 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <gtest/gtest.h>
#include <math.h>
extern "C" {
#include <multirow/double.h>
#include <multirow/util.h>
}
TEST(DoubleTest, double_cmp_test)
{
EXPECT_EQ(0, DOUBLE_cmp(1/3.0, 2/6.0));
EXPECT_EQ(1, DOUBLE_cmp(5.0, 3.0));
EXPECT_EQ(0, DOUBLE_cmp(5.0, 5.0));
EXPECT_EQ(-1, DOUBLE_cmp(5.0, 10.0));
EXPECT_EQ(1, DOUBLE_cmp(5.1, 5.0));
EXPECT_EQ(0, DOUBLE_cmp(1000000005, 1000000003));
EXPECT_EQ(0, DOUBLE_cmp(1000000005, 1000000005));
EXPECT_EQ(0, DOUBLE_cmp(1000000005, 1000000007));
EXPECT_EQ(1, DOUBLE_cmp(0.00005, 0.00003));
EXPECT_EQ(0, DOUBLE_cmp(0.00005, 0.00005));
EXPECT_EQ(-1, DOUBLE_cmp(0.00005, 0.00007));
EXPECT_EQ(0, DOUBLE_cmp(0.5 * EPSILON, 0.3 * EPSILON));
EXPECT_EQ(0, DOUBLE_cmp(0.5 * EPSILON, 0.5 * EPSILON));
EXPECT_EQ(0, DOUBLE_cmp(0.5 * EPSILON, 0.7 * EPSILON));
EXPECT_EQ(0, DOUBLE_cmp(0.0, 0.0));
EXPECT_EQ(0, DOUBLE_cmp(0.0, -0.0));
EXPECT_EQ(0, DOUBLE_cmp(-0.0, -0.0));
EXPECT_EQ(0, DOUBLE_cmp(0.1 * EPSILON, 0.0));
EXPECT_EQ(0, DOUBLE_cmp(0.0, 0.1 * EPSILON));
EXPECT_EQ(0, DOUBLE_cmp(-0.1 * EPSILON, 0.0));
EXPECT_EQ(0, DOUBLE_cmp(0.0, -0.1 * EPSILON));
EXPECT_EQ(0, DOUBLE_cmp(DBL_MAX, DBL_MAX));
EXPECT_EQ(1, DOUBLE_cmp(DBL_MAX, -DBL_MAX));
EXPECT_EQ(-1, DOUBLE_cmp(-DBL_MAX, DBL_MAX));
EXPECT_EQ(1, DOUBLE_cmp(DBL_MAX, DBL_MAX / 2));
EXPECT_EQ(1, DOUBLE_cmp(DBL_MAX, -DBL_MAX / 2));
EXPECT_EQ(-1, DOUBLE_cmp(-DBL_MAX, DBL_MAX / 2));
EXPECT_EQ(0, DOUBLE_cmp(INFINITY, INFINITY));
EXPECT_EQ(1, DOUBLE_cmp(INFINITY, -INFINITY));
EXPECT_EQ(-1, DOUBLE_cmp(-INFINITY, INFINITY));
EXPECT_EQ(0, DOUBLE_cmp(-INFINITY, -INFINITY));
EXPECT_EQ(1, DOUBLE_cmp(INFINITY, DBL_MAX));
EXPECT_EQ(-1, DOUBLE_cmp(-INFINITY, DBL_MIN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, 0.0f));
EXPECT_EQ(1, DOUBLE_cmp(-0.0f, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, -0.0f));
EXPECT_EQ(1, DOUBLE_cmp(0.0f, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, INFINITY));
EXPECT_EQ(1, DOUBLE_cmp(INFINITY, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, -INFINITY));
EXPECT_EQ(1, DOUBLE_cmp(-INFINITY, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, DBL_MAX));
EXPECT_EQ(1, DOUBLE_cmp(DBL_MAX, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, -DBL_MAX));
EXPECT_EQ(1, DOUBLE_cmp(-DBL_MAX, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, DBL_MIN));
EXPECT_EQ(1, DOUBLE_cmp(DBL_MIN, NAN));
EXPECT_EQ(1, DOUBLE_cmp(NAN, -DBL_MIN));
EXPECT_EQ(1, DOUBLE_cmp(-DBL_MIN, NAN));
}
TEST(DoubleTest, double_iszero_test)
{
EXPECT_TRUE(DOUBLE_iszero(0.0));
EXPECT_TRUE(DOUBLE_iszero(-0.0));
EXPECT_TRUE(DOUBLE_iszero(0.1 * EPSILON));
EXPECT_TRUE(DOUBLE_iszero(-0.1 * EPSILON));
}
TEST(DoubleTest, double_eq_test)
{
EXPECT_TRUE(DOUBLE_eq(1.0, 1.0));
EXPECT_FALSE(DOUBLE_eq(1.0, 3.0));
EXPECT_TRUE(DOUBLE_eq(1/3.0, 2/6.0));
}
TEST(DoubleTest, double_max_test)
{
EXPECT_EQ(0.1, DOUBLE_max(0.0, 0.1));
EXPECT_EQ(0.5, DOUBLE_max(0.5, 0.1));
}
TEST(DoubleTest, double_geq_test)
{
EXPECT_TRUE(DOUBLE_geq(1.0, 0.5));
EXPECT_TRUE(DOUBLE_geq(0.46000000, 0.45953916));
}
TEST(DoubleTest, double_to_rational_test)
{
int rval = 0;
Rational r;
rval = DOUBLE_to_rational(M_PI, 1, r);
abort_if(rval, "DOUBLE_to failed");
EXPECT_EQ(3, r->num);
EXPECT_EQ(1, r->den);
rval = DOUBLE_to_rational(M_PI, 10, r);
abort_if(rval, "DOUBLE_to failed");
EXPECT_EQ(22, r->num);
EXPECT_EQ(7, r->den);
rval = DOUBLE_to_rational(M_PI, 100, r);
abort_if(rval, "DOUBLE_to failed");
EXPECT_EQ(311, r->num);
EXPECT_EQ(99, r->den);
rval = DOUBLE_to_rational(M_PI, 1000, r);
abort_if(rval, "DOUBLE_to failed");
EXPECT_EQ(355, r->num);
EXPECT_EQ(113, r->den);
rval = DOUBLE_to_rational(M_PI, 10000, r);
abort_if(rval, "DOUBLE_to failed");
EXPECT_EQ(355, r->num);
EXPECT_EQ(113, r->den);
rval = DOUBLE_to_rational(-M_PI, 10000, r);
abort_if(rval, "DOUBLE_to failed");
EXPECT_EQ(-355, r->num);
EXPECT_EQ(113, r->den);
CLEANUP:
if(rval) FAIL();
}
TEST(DoubleTest, gcd_test)
{
EXPECT_EQ(1, gcd(11, 7));
EXPECT_EQ(2, gcd(6, 4));
EXPECT_EQ(2, gcd(4, 6));
EXPECT_EQ(5, gcd(5, 10));
}

664
multirow/tests/geometry-test.cpp
Unescape View File

@ -0,0 +1,664 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <gtest/gtest.h>
extern "C" {
#include <math.h>
#include <multirow/util.h>
#include <multirow/geometry.h>
}
TEST(GeometryTest, interpolate_coefficients_1)
{
int rval = 0;
double r1[] = {1.0, 2.0};
double r2[] = {3.0, 1.0};
double p[] = {2.0, 2.0};
double lambda;
rval = scale_vector_to_line(r1, r2, p, &lambda);
abort_if(rval, "scale_vector_to_line failed");
EXPECT_DOUBLE_EQ(5 / 6.0, lambda);
CLEANUP:
if (rval) FAIL();
}
TEST(GeometryTest, interpolate_coefficients_2)
{
int rval = 0;
double r1[] = {1.0, 3.0};
double r2[] = {1e90, 1e90};
double p1[] = {2.0, 3.0};
double p2[] = {1.0, 1.0};
double p3[] = {1.0, 3.0};
double lambda1;
double lambda2;
double lambda3;
rval = scale_vector_to_line(r1, r2, p1, &lambda1);
abort_if(rval, "scale_vector_to_line failed");
rval = scale_vector_to_line(r1, r2, p2, &lambda2);
abort_if(rval, "scale_vector_to_line failed");
rval = scale_vector_to_line(r1, r2, p3, &lambda3);
abort_if(rval, "scale_vector_to_line failed");
EXPECT_DOUBLE_EQ(2.0, lambda1);
EXPECT_DOUBLE_EQ(INFINITY, lambda2);
EXPECT_DOUBLE_EQ(1.0, lambda3);
rval = scale_vector_to_line(r2, r1, p1, &lambda1);
abort_if(rval, "scale_vector_to_line failed");
rval = scale_vector_to_line(r2, r1, p2, &lambda2);
abort_if(rval, "scale_vector_to_line failed");
rval = scale_vector_to_line(r2, r1, p3, &lambda3);
abort_if(rval, "scale_vector_to_line failed");
EXPECT_DOUBLE_EQ(2.0, lambda1);
EXPECT_DOUBLE_EQ(INFINITY, lambda2);
EXPECT_DOUBLE_EQ(1.0, lambda3);
CLEANUP:
if (rval) FAIL();
}
TEST(GeometryTest, interpolate_coefficients_3)
{
int rval = 0;
double r1[] = {-1.0, 1.0};
double r2[] = {-1.0, -1.0};
double p1[] = {-1.0, 0.0};
double lambda1;
rval = scale_vector_to_line(r1, r2, p1, &lambda1);
abort_if(rval, "scale_vector_to_line failed");
EXPECT_DOUBLE_EQ(1.0, lambda1);
CLEANUP:
if (rval) FAIL();
}
TEST(GeometryTest, chull_2d_test)
{
int rval = 0;
int npoints = 9;
double points[] =
{
0.0, 0.0,
0.0, 1.0,
0.0, 2.0,
1.0, 0.0,
1.0, 1.0,
1.0, 2.0,
2.0, 0.0,
2.0, 1.0,
2.0, 2.0
};
double vertices[100];
int nvertices;
rval = chull_2d(points, npoints, vertices, &nvertices);
abort_if(rval, "chull_2d failed");
for(int i=0; i<nvertices; i++)
log_debug(" %.2lf %.2lf\n", vertices[2*i], vertices[2*i+1]);
CLEANUP:
if(rval) FAIL();
}
//TEST(GeometryTest, sort_rays_test)
//{
// double rays[] = {1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0};
// int var_to_ray[] = {0, 0, 1, 2, 3};
// double bounds[] = {1.0, 2.0, 3.0, 4.0};
//
// sort_rays(rays, 4, bounds, var_to_ray, 5);
//
// EXPECT_EQ(rays[0], -1.0);
// EXPECT_EQ(rays[1], -1.0);
//
// EXPECT_EQ(rays[2], -1.0);
// EXPECT_EQ(rays[3], 1.0);
//
// EXPECT_EQ(rays[4], 1.0);
// EXPECT_EQ(rays[5], 1.0);
//
// EXPECT_EQ(rays[6], 1.0);
// EXPECT_EQ(rays[7], -1.0);
//
// EXPECT_EQ(var_to_ray[0], 2);
// EXPECT_EQ(var_to_ray[1], 2);
// EXPECT_EQ(var_to_ray[2], 1);
// EXPECT_EQ(var_to_ray[3], 3);
// EXPECT_EQ(var_to_ray[4], 0);
//
// EXPECT_EQ(bounds[0], 4.0);
// EXPECT_EQ(bounds[1], 2.0);
// EXPECT_EQ(bounds[2], 1.0);
// EXPECT_EQ(bounds[3], 3.0);
//}
//
//TEST(GeometryTest, scale_rays)
//{
// int rval = 0;
//
// double rays[] = {1.0, 1.0, 2.0, 1.0};
// int nrays = 2;
// double scale[] = {5.0, 2.0};
//
// rval = scale_rays(rays, nrays, scale);
// abort_if(rval, "scale_rays failed");
//
// EXPECT_DOUBLE_EQ(rays[0], 5.0);
// EXPECT_DOUBLE_EQ(rays[1], 5.0);
// EXPECT_DOUBLE_EQ(rays[2], 4.0);
// EXPECT_DOUBLE_EQ(rays[3], 2.0);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, scale_cone_to_point_test)
//{
// int rval = 0;
//
// double r1[2] = { 1.0, 0.0 };
// double r2[2] = { 0.0, 1.0 };
//
// double p1[2] = { 1.0, 1.0 };
// double p2[2] = { 1.0, 0.0 };
// double p3[2] = { -1.0, 1.0 };
// double p4[2] = { -1.0, 0.0 };
//
// double lambda;
//
// rval = scale_cone_to_point(r1, r2, p1, &lambda);
// abort_if(rval, "scale_cone_to_point failed");
// EXPECT_DOUBLE_EQ(2.0, lambda);
//
// rval = scale_cone_to_point(r1, r2, p2, &lambda);
// abort_if(rval, "scale_cone_to_point failed");
// EXPECT_DOUBLE_EQ(1.0, lambda);
//
// rval = scale_cone_to_point(r1, r2, p3, &lambda);
// abort_if(rval, "scale_cone_to_point failed");
// EXPECT_DOUBLE_EQ(0.0, lambda);
//
// rval = scale_cone_to_point(r1, r2, p4, &lambda);
// abort_if(rval, "scale_cone_to_point failed");
// EXPECT_DOUBLE_EQ(-1.0, lambda);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, shear_cone_to_point_test)
//{
// int rval = 0;
//
// double r1[2] = { 0.0, 1.0 };
// double r2[2] = { 1.0, 0.0 };
//
// double p1[2] = { 1.0, 0.5 };
// double p2[2] = { 1.5, 0.5 };
// double p3[2] = { 0.5, 0.0 };
// double p4[2] = { 1.0, 1.0 };
// double p5[2] = { 1.0, 2.0 };
//
// double lambda;
//
// rval = shear_cone_to_point(r1, r2, p1, &lambda);
// abort_if(rval, "shear_cone_to_point failed");
// EXPECT_DOUBLE_EQ(2.0, lambda);
//
// rval = shear_cone_to_point(r1, r2, p2, &lambda);
// abort_if(rval, "shear_cone_to_point failed");
// EXPECT_DOUBLE_EQ(3.0, lambda);
//
// rval = shear_cone_to_point(r1, r2, p3, &lambda);
// abort_if(rval, "shear_cone_to_point failed");
// EXPECT_DOUBLE_EQ(0.5, lambda);
//
// rval = shear_cone_to_point(r1, r2, p4, &lambda);
// abort_if(rval, "shear_cone_to_point failed");
// EXPECT_DOUBLE_EQ(INFINITY, lambda);
//
// rval = shear_cone_to_point(r1, r2, p5, &lambda);
// abort_if(rval, "shear_cone_to_point failed");
// EXPECT_DOUBLE_EQ(-1.0, lambda);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bound_for_point_test_1)
//{
// int rval = 0;
//
// int nrays = 4;
// double f[] = { 0.75, 0.75 };
// double rays[] = { 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, -1.0, 0.0 };
// double bounds[] = { INFINITY, INFINITY, INFINITY, INFINITY };
//
// double p1[] = { 0.0, 0.0 };
// double p2[] = { 0.0, 1.0 };
// double p3[] = { 1.0, 0.0 };
// double p4[] = { 1.0, 1.0 };
//
// double epsilon;
// double v1[2], v2[2];
// int i1, i2;
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p1, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(1.5, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p2, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(1.0, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p3, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(1.0, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p4, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(0.5, epsilon);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bound_for_point_test_2)
//{
// int rval = 0;
//
// int nrays = 4;
//
// double f[] = { 0.75, 0.75 };
// double rays[] = { 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, -1.0, 0.0 };
// double bounds[] = { 0.5, 0.5, INFINITY, INFINITY };
//
// double p1[] = { 0.0, 0.0 };
// double p2[] = { 0.0, 1.0 };
// double p3[] = { 1.0, 0.0 };
// double p4[] = { 1.0, 1.0 };
//
// double epsilon;
// double v1[2], v2[2];
// int i1, i2;
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p1, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(1.5, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p2, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(1.5, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p3, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(1.5, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p4, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(INFINITY, epsilon);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bound_for_point_test_3)
//{
// int rval = 0;
//
// int nrays = 4;
//
// double f[] = { 0.75, 0.75 };
// double rays[] = { 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, -1.0, 0.0 };
// double bounds[] = { 0.5, 0.5, 1.5, 1.5 };
//
// double p1[] = { 0.0, 0.0 };
// double p2[] = { 0.0, 1.0 };
// double p3[] = { 1.0, 0.0 };
// double p4[] = { 1.0, 1.0 };
//
// double epsilon;
// double v1[2], v2[2];
// int i1, i2;
//
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p1, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(INFINITY, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p2, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(INFINITY, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p3, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(INFINITY, epsilon);
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p4, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(INFINITY, epsilon);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bound_for_point_test_4)
//{
// int rval = 0;
//
// int nrays = 4;
//
// double f[] = { 0.25, 0.5 };
// double rays[] = { 1.0, 1.0, 0.0, -1.0, -1.0, 1.0 };
// double bounds[] = { 0.5, INFINITY, 0.5, INFINITY };
//
// double p1[] = { 0.0, 0.0 };
//
// double epsilon;
// double v1[2], v2[2];
// int i1, i2;
//
// rval = compute_bound_for_point(rays, bounds, nrays, f, p1, &epsilon, v1, v2, &i1, &i2);
// abort_if(rval, "compute_bound_for_point failed");
// EXPECT_DOUBLE_EQ(1.5, epsilon);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bounds_2d_test_1)
//{
// int rval = 0;
//
// int nrays = 4;
// double f[] = { 0.75, 0.75 };
// double rays[] = { 0.0, 1.0, 1.0, 0.0, 0.0, -1.0, -1.0, 0.0 };
//
// double bounds[4];
//
// rval = compute_bounds_2d(rays, nrays, f, bounds);
// abort_if(rval, "compute_bounds_2d failed");
//
// EXPECT_DOUBLE_EQ(0.5, bounds[0]);
// EXPECT_DOUBLE_EQ(0.5, bounds[1]);
// EXPECT_DOUBLE_EQ(1.5, bounds[2]);
// EXPECT_DOUBLE_EQ(1.5, bounds[3]);
//
//CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bounds_2d_test_2)
//{
// int rval = 0;
//
// int nrays = 3;
// double f[] = { 0.25, 0.5 };
// double rays[] = { 1.0, 1.0, 0.0, -1.0, -1.0, 1.0 };
//
// double bounds[3];
//
// rval = compute_bounds_2d(rays, nrays, f, bounds);
// abort_if(rval, "compute_bounds_2d failed");
//
// EXPECT_DOUBLE_EQ(0.5, bounds[0]);
// EXPECT_DOUBLE_EQ(1.5, bounds[1]);
// EXPECT_DOUBLE_EQ(0.5, bounds[2]);
//
//CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bounds_2d_test_3)
//{
// int rval = 0;
//
// int nrays = 5;
// double f[] = {1 / 4.0, 3 / 4.0};
// double rays[] = {-2 / 5.0, 5 / 7.0, 0.0, 1.0, 1.0, 1.0, 4 / 5.0, -2 / 3.0,
// -1.0, 0.0};
//
// double bounds[5];
//
// rval = compute_bounds_2d(rays, nrays, f, bounds);
// abort_if(rval, "compute_bounds_2d failed");
//
// EXPECT_DOUBLE_EQ(23 / 50.0, bounds[0]);
// EXPECT_DOUBLE_EQ(23 / 42.0, bounds[1]);
// EXPECT_DOUBLE_EQ(9 / 11.0, bounds[2]);
// EXPECT_DOUBLE_EQ(9 / 11.0, bounds[3]);
// EXPECT_DOUBLE_EQ(23 / 50.0, bounds[4]);
//
// CLEANUP:
// if (rval) FAIL();
//}
//
//
//TEST(GeometryTest, compute_bounds_2d_test_4)
//{
// int rval = 0;
//
// int nrays = 5;
// double f[] = {1 / 2.0, 1 / 2.0};
// double rays[] = {-1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, -1.0, -1.0, -1.0 };
//
// double bounds[5];
//
// rval = compute_bounds_2d(rays, nrays, f, bounds);
// abort_if(rval, "compute_bounds_2d failed");
//
// EXPECT_DOUBLE_EQ(0.5, bounds[0]);
// EXPECT_DOUBLE_EQ(0.5, bounds[1]);
// EXPECT_EQ(GREEDY_BIG_E, bounds[2]);
// EXPECT_DOUBLE_EQ(0.5, bounds[3]);
// EXPECT_DOUBLE_EQ(0.5, bounds[4]);
//
// CLEANUP:
// if (rval) FAIL();
//}
//
//TEST(GeometryTest, compute_bounds_2d_test_5)
//{
// int rval = 0;
//
// int nrays = 6;
// double f[] = {1 / 2.0, 1 / 2.0};
// double rays[] = {-1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, -1.0, -1.0, -1.0 };
//
// double bounds[6];
//
// rval = compute_bounds_2d(rays, nrays, f, bounds);
// abort_if(rval, "compute_bounds_2d failed");
//
// EXPECT_DOUBLE_EQ(0.5, bounds[0]);
// EXPECT_DOUBLE_EQ(1.0, bounds[1]);
// EXPECT_DOUBLE_EQ(0.5, bounds[2]);
// EXPECT_DOUBLE_EQ(1.0, bounds[3]);
// EXPECT_DOUBLE_EQ(0.5, bounds[4]);
// EXPECT_DOUBLE_EQ(0.5, bounds[5]);
//
// CLEANUP:
// if (rval) FAIL();
//}
//
////TEST(GeometryTest, next_lattice_point_test)
////{
//// int rval = 0;
////
//// LatticeSequence seq;
//// lattice_sequence_init(&seq);
////
//// double lb[] = { -9.0, -1.0 };
//// double ub[] = { 3.0, 4.0 };
////
////
//// while(!seq.eol)
//// {
//// printf("%3d %3d\n", seq.i, seq.j);
//// next_lattice_point(&seq, lb, ub);
//// }
////
//// CLEANUP:
//// if(rval) FAIL();
////}
//
//
//TEST(GeometryTest, get_plane_test)
//{
// int rval = 0;
//
// double v0[] = { -1.0, 1.0 };
// double v1[] = { 1.0, 1.0 };
// double v2[] = { 0.0, 1.0 };
// double v3[] = { 1.0, 0.0 };
//
// double plane[2];
//
// rval = get_plane(v0, v1, plane);
// abort_if(rval, "get_plane failed");
// EXPECT_DOUBLE_EQ(0.0, plane[0]);
// EXPECT_DOUBLE_EQ(1.0, plane[1]);
//
// rval = get_plane(v1, v0, plane);
// abort_if(rval, "get_plane failed");
// EXPECT_DOUBLE_EQ(0.0, plane[0]);
// EXPECT_DOUBLE_EQ(1.0, plane[1]);
//
// rval = get_plane(v2, v3, plane);
// abort_if(rval, "get_plane failed");
// EXPECT_DOUBLE_EQ(1.0, plane[0]);
// EXPECT_DOUBLE_EQ(1.0, plane[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, generate_split_test)
//{
// int rval = 0;
//
// double f[2] = { 0.5, 0.0 };
// double d[2] = { 1.0, 1.0 };
//
// double pi[2], pi_zero;
//
// rval = generate_split(f, d, pi, &pi_zero, 1000);
// abort_if(rval, "generate_split failed");
//
// EXPECT_DOUBLE_EQ(0.0, pi_zero);
// EXPECT_DOUBLE_EQ(1.0, pi[0]);
// EXPECT_DOUBLE_EQ(-1.0, pi[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, generate_split_test_2)
//{
// int rval = 0;
//
// double f[2] = { 0.5, 0.25 };
// double d[2] = { 9/15.0, 13/21.0 };
//
// double pi[2], pi_zero;
//
// rval = generate_split(f, d, pi, &pi_zero, 1000);
// abort_if(rval, "generate_split failed");
//
// EXPECT_DOUBLE_EQ(16.0, pi_zero);
// EXPECT_DOUBLE_EQ(65.0, pi[0]);
// EXPECT_DOUBLE_EQ(-63.0, pi[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, generate_split_test_3)
//{
// int rval = 0;
//
// double f[2] = { 0.0, 0.5 };
// double d[2] = { -1.0, 0.0 };
//
// double pi[2], pi_zero;
//
// rval = generate_split(f, d, pi, &pi_zero, 1000);
// abort_if(rval, "generate_split failed");
//
// EXPECT_DOUBLE_EQ(0.0, pi_zero);
// EXPECT_DOUBLE_EQ(0.0, pi[0]);
// EXPECT_DOUBLE_EQ(1.0, pi[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//
//TEST(GeometryTest, generate_split_test_4)
//{
// int rval = 0;
//
// double f[2] = { 0.5, 0.5 };
// double d[2] = { 0.0, 1.70 };
//
// double pi[2], pi_zero;
//
// rval = generate_split(f, d, pi, &pi_zero, 1000);
// abort_if(rval, "generate_split failed");
//
// EXPECT_DOUBLE_EQ(0.0, pi_zero);
// EXPECT_DOUBLE_EQ(1.0, pi[0]);
// EXPECT_DOUBLE_EQ(0.0, pi[1]);
//
// CLEANUP:
// if(rval) FAIL();
//}
//TEST(GeometryTest, find_linear_combination)
//{
// double p[] = {1.0, 1.0};
// double r1[] = {1 / 3.0, 1 / 5.0};
// double r2[] = {-1.0, 2 / 3.0};
// double lambda1, lambda2;
//
// find_linear_combination(r1, r2, p, &lambda1, &lambda2);
//
// EXPECT_NEAR(3.94736842105263, lambda1, 1e-12);
// EXPECT_NEAR(0.315789473684211, lambda2, 1e-12);
//}