diff --git a/Code/hull.cpp b/Code/hull.cpp
index 0198facce2a64dbc73e0c3164a806bbab19af69f..b0f9b693e17fa3eab4c3f538f996df919e94ea0e 100644
--- a/Code/hull.cpp
+++ b/Code/hull.cpp
@@ -38,11 +38,6 @@ struct vec {
float length() const {
return sqrt(length2()) ;
}
-
- bool lexicoinf(const vec& rhs) const {
- if(x == rhs.x) return y < rhs.y ;
- return x < rhs.x ;
- }
} ;
//}}}
@@ -53,29 +48,10 @@ int sign(float x) {
return 0 ;
}
-int orient(const vec& v0_in, const vec& v1_in) {
- vec v0, v1 ;
- int factor = 1 ;
- if(v0_in.lexicoinf(v1_in)) {
- v0 = v0_in ;
- v1 = v1_in ;
- } else {
- v0 = v1_in ;
- v1 = v0_in ;
- factor = -1 ;
- }
-
+int orient(const vec& v0, const vec& v1) {
//determinant
float d = v0.x*v1.y - v0.y*v1.x ;
- return d ;
-
- if(d != 0) return factor*sign(d) ;
-
- //perturbation
- if(v1.y != 0) return factor*sign(v1.y) ;
- float a = v0.x - v0.y - v1.x ;
- if (a!= 0) return factor*sign(a) ;
- return factor*1 ;
+ return sign(d) ;
}
struct vec_compare {
@@ -94,8 +70,10 @@ struct vec_compare {
size_t leftmost(const std::vector<vec>& points) {
size_t min = 0 ;
+ //iterato on all the points
for(size_t i = 1; i < points.size(); ++i) {
if(points[i].x < points[min].x) {
+ //update the leftmost
min = i ;
}
}
@@ -139,21 +117,27 @@ void svg_hull(
float salt = 0
)
{
+ //open file for export
std::ofstream file(filename) ;
SvgPad svg(file, 1024, 1024) ;
svg.open() ;
+ //plot data points
for(const vec& p : points) {
svg.solid_point(p.x, p.y) ;
}
- size_t hs = hull.size() ;
+ //randomness in case salt is added for hull debugging
std::mt19937 mt ;
std::uniform_real_distribution<float> rf(-salt, salt) ;
+
+ //plot hull
+ size_t hs = hull.size() ;
for(size_t i = 0; i < hs; ++i) {
- vec salti(0,0) ;
- vec saltip(0,0) ;
+ //salt to perturb the line endpoints and avoid overlaps
+ vec salti(0,0) ;
+ vec saltip(0,0) ;
if(salt > 0) {
salti.x = rf(mt) ;
salti.y = rf(mt) ;
@@ -161,7 +145,10 @@ void svg_hull(
saltip.y = rf(mt) ;
}
+ //hull points
svg.contour_point(hull[i].x, hull[i].y) ;
+
+ //hull sides with salt
svg.line(
hull[i].x + salti.x,
hull[i].y + salti.y,
@@ -170,23 +157,31 @@ void svg_hull(
) ;
}
+ //close file
svg.close() ;
}
+//utility to illustrate angle sorting
void svg_sort(
const std::string& filename,
const std::vector<vec>& points)
{
+ //open svg file for export
std::ofstream file(filename) ;
SvgPad svg(file, 1024, 1024) ;
svg.open() ;
+ //colors for the linear gradient of lines
float col0[3] = { 1., 0.9, 0.3 } ;
float col1[3] = { 0.35, 0.2, 0.4 } ;
+
+ //iterate on every point except the leftmost at 0
for(size_t i = 1; i < points.size(); ++i) {
+ //hue
float b = ((float) i-1) / (points.size() - 1) ;
float a = 1-b ;
svg.stroke(a*col0[0] + b*col1[0], a*col0[1] + b*col1[1], a*col0[2] + b*col1[2]) ;
+ //line defining the sorting angle
svg.line(points[0].x, points[0].y, points[i].x, points[i].y) ;
}
@@ -236,10 +231,9 @@ int main() {
aligned_points.emplace_back(0, 0.8) ;
aligned_points.emplace_back(0, 0.5) ;
+ //test every permutation of the input
int indices[4] = {0, 1, 2, 3} ;
int permutation_index = 0 ;
-
- //test every permutation of the input
do {
std::vector<vec> points ;
for(unsigned int i = 0; i < 4; ++i) {
@@ -252,7 +246,7 @@ int main() {
std::vector<vec> hull ;
compute_hull(points, hull) ;
- //export
+ //export with salt for readability
std::stringstream ss ;
ss << "/tmp/hull_aligned_" << permutation_index << ".svg" ;
svg_hull(ss.str(), points, hull, 0.02) ;
diff --git a/Code/precision.cpp b/Code/precision.cpp
index f9bca1fb07470c2842ddb59d42dc3156851d1824..22732da4de4c7d0ad72795b3fcf6419555dbc4de 100644
--- a/Code/precision.cpp
+++ b/Code/precision.cpp
@@ -16,6 +16,7 @@ using VecIX = Eigen::Matrix<interval, Eigen::Dynamic, 1> ;
using MatRX = Eigen::Matrix<rational, Eigen::Dynamic, Eigen::Dynamic> ;
using VecRX = Eigen::Matrix<rational, Eigen::Dynamic, 1> ;
+//fix problems when using intervals in Eigen
namespace Eigen {
namespace internal {
template<typename X, typename S, typename P>
@@ -30,62 +31,80 @@ namespace Eigen {
}
}
-//command line
-char* get_cmd_option( char ** begin, char ** end,
- const std::string & option
- ) {
- char ** itr = std::find(begin, end, "--" + option);
- if (itr != end && ++itr != end)
- {
- return *itr;
- }
- return 0;
+//{{{ command line
+
+//get the token immediatly after an option given by name
+char* get_cmd_option(
+ char ** begin, char ** end,
+ const std::string & option
+ )
+{
+ //look for the desired option
+ char ** itr = std::find(begin, end, "--" + option);
+ if (itr != end && ++itr != end)
+ {
+ //the option was found
+ return *itr;
+ }
+ return 0;
}
-bool cmd_option_exists( char** begin, char** end,
- const std::string& option
- ) {
- return std::find(begin, end, "--" + option) != end;
+//check the existence of an option
+bool cmd_option_exists(
+ char** begin, char** end,
+ const std::string& option
+ )
+{
+ //look for the desired option
+ return std::find(begin, end, "--" + option) != end;
}
+//try to fill a provided variable with a command line option
template<typename T>
bool get_cmd_value(
char ** begin,
char** end,
const std::string& option,
- T* output
- ) {
+ T& output
+ )
+{
if(cmd_option_exists(begin, end, option)) {
+ //the option exists, get the next token
std::stringstream ss ;
ss << get_cmd_option(begin, end, option) ;
+ //try to parse the token as the provided type
T t ;
ss >> t ;
if(!ss.fail()) {
- *output = t ;
+ //the parsing is allowed, fill the output
+ output = t ;
return true ;
}
}
+ //no parsing was possible, leave the variable at its initial value
return false ;
}
+//}}}
+
int main(int argc, char** argv) {
//test variables
- unsigned int dim = 10 ;
- get_cmd_value(argv, argv + argc, "dim", &dim) ;
+ unsigned int dim = 20 ;
+ get_cmd_value(argv, argv + argc, "dim", dim) ;
//alea
std::random_device seed ;
std::mt19937_64 mt(seed()) ;
std::uniform_real_distribution<double> random_double(-1, 1) ;
- //initialization
+ //double matrix
MatDX matd(dim, dim) ;
- VecDX vecd(dim) ;
+ //interval matrix
MatIX mati(dim, dim) ;
- VecIX veci(dim) ;
+ //rational matrix
MatRX matr(dim, dim) ;
- VecRX vecr(dim) ;
+ //randomly fill the matrices with the same values
for(unsigned int i = 0; i < dim; ++i) {
for(unsigned int j = 0; j < dim; ++j) {
double d = random_double(mt) ;
@@ -93,40 +112,52 @@ int main(int argc, char** argv) {
mati(i, j) = d ;
matr(i, j) = d ;
}
- double d = random_double(mt) ;
- vecd(i) = d ;
- veci(i) = d ;
- vecr(i) = d ;
}
+ //compute the exact determinant in rational
rational detr = matr.fullPivLu().determinant() ;
std::cout << "==========" << std::endl ;
{
+ //testing fullPivLu, supposedly slower but more precise
std::cout << "solve fullPivLu" << std::endl ;
+ //double determinant
double detd = matd.fullPivLu().determinant() ;
try {
+ //interval determinant
interval deti = mati.fullPivLu().determinant() ;
+ //computed double determinant
std::cout << detd
+ //interval bounds
<< " in [" << deti.lower() << " -- " << deti.upper()
- << "] error " << detd - (double) detr
+ //error wrt. ground truth
+ << "] error " << (double) (detd - detr)
+ //width of the error interval
<< " wrt. " << deti.upper() - deti.lower() << std::endl ;
} catch(std::exception& e) {
+ //using interval yielded an exception, fallback
std::cout << "interval computation failed" << std::endl ;
std::cout << detd << " error " << detd - (double) detr << std::endl ;
}
std::cout << "==========" << std::endl ;
}
{
+ //testing fullPivLu, supposedly faster but less precise
std::cout << "solve partialPivLu" << std::endl ;
double detd = matd.partialPivLu().determinant() ;
try {
+ //interval determinant
interval deti = mati.partialPivLu().determinant() ;
+ //computed double determinant
std::cout << detd
+ //interval bounds
<< " in [" << deti.lower() << " -- " << deti.upper()
+ //error wrt. ground truth
<< "] error " << detd - (double) detr
+ //width of the error interval
<< " wrt. " << deti.upper() - deti.lower() << std::endl ;
} catch(std::exception& e) {
+ //using interval yielded an exception, fallback
std::cout << "interval computation failed" << std::endl ;
std::cout << detd << " error " << detd - (double) detr << std::endl ;
}