Skip to content
Snippets Groups Projects
Commit 5e62402e authored by Vincent Nivoliers's avatar Vincent Nivoliers
Browse files

more coments

parent 8b0cf459
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Code/hull.cpp
+ 26
− 32
View file @ 5e62402e
......@@ -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) ;
......
Code/precision.cpp
+ 59
− 28
View file @ 5e62402e
......@@ -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 ;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment