phynexis-deps/cgal
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Files
62908b7fafdcbb79175bb00b0281fa4d650ca874
cgal/Kinetic_data_structures/test/Polynomial_kernel/include/Check_solver.h
T
Eric Berberich 681f282d77 Fix by Daniel Russel to make testsuite greener 
Raise "NUMERICAL ISSUE" instead of "WARNING". The latter made some 
rows of the testsuite completely yellow. Fix has been tested locally.
2012-10-09 12:55:09 +00:00

473 lines
13 KiB
C++
Raw Blame History

#ifndef CHECK_SOLVER_H
#define CHECK_SOLVER_H
#include <CGAL/Polynomial/internal/polynomial_generators.h>
#include <CGAL/Timer.h>
#include <vector>
#include <iterator>
struct True{};
struct False{};
template <class K, class CE=False, bool CA=true>
class Check_solver
{
public:
typedef typename K::Root_stack Root_stack;
typedef typename K::Root_container Root_container;
typedef std::vector<double> DV;
typedef typename K::Root Rt;
typedef typename K::Function Fn;
typedef typename Fn::NT NT;
typedef K Kernel;
K k_;
bool verbose;
typename K::Construct_function cf_;
double total_time;
unsigned int num_roots;
double ce_time;
unsigned int nce;
Check_solver(const K& k, bool verbose): k_(k), verbose(verbose), cf_(k_.construct_function_object()){}
/*void time_estimate(const Fn &, const Rt& , False){
}
void time_estimate(const Fn &q, const Rt& start, True){
Root_stack s= k_.root_stack_object(q, start, std::numeric_limits<Rt>::infinity());
CGAL::Timer timer;
timer.start();
do {
s.compute_estimate(start);
++nce;
} while (timer.time() <3 && !verbose);
timer.stop();
ce_time= timer.time();
}*/
template <class It>
void check_polynomial(const Fn &q, It roots_b, It roots_e,
const Rt& start = Rt(-1000.0),
const Rt& end = Rt(1000.0) ) {
CGAL::Timer timer;
//typename K::Is_even_multiplicity iem = k_.is_even_multiplicity_object(q);
//int current_root=0;
if (verbose) {
std::cout << "Polynomial: " << q << std::endl;
std::cout << "Window: [";
std::cout<< CGAL::to_double(start);
std::cout << ", ";
std::cout << CGAL::to_double(end);
std::cout << ")" << std::endl;
std::cout << "Solution: ";
for (It c= roots_b; c!= roots_e; ++c) {
/*if (CGAL::to_double(start) >= a[i]-.00001){
++current_root;
} else if (CGAL::to_double(end) <= a[i]){
} else {*/
std::cout << "<" << *c << "> ";
//}
}
std::cout << std::endl;
std::cout << "Solver: ";
}
std::vector<Rt> roots;
DV comp;
int reps=0;
timer.start();
int total_roots=0;
do {
Root_stack s= k_.root_stack_object(q, start, end);
Rt last_root= -Rt(1000.0);
while (!s.empty()) {
if (reps==0 && verbose) {
//comp.push_back(CGAL::to_double(r));
if (last_root > s.top()) {
std::cerr << "NUMERIC ISSUE last root was " << last_root << " and current root is "
<< s.top() << std::endl;
}
CGAL_assertion(last_root<= s.top());
last_root=s.top();
//Rt cur= s.top();
//if (cur != std::numeric_limits<Rt>::infinity()) {
/*if (cur == std::numeric_limits<Rt>::infinity() || cur > 10000){
int i=0;
}*/
roots.push_back(s.top());
std::cout << "<" << CGAL::to_double(s.top());
//if (iem(s.top())) { std::cout <<"E";}
std::cout << "> " << std::flush;
/*}
else {
}*/
}
s.pop();
++total_roots;
}
++total_roots;
if (reps==0 && verbose && CA) {
std::vector<bool> taken_maple(roots_e-roots_b, false), taken_solver(roots.size(), false);
for (unsigned int i=0; i< roots.size(); ++i) {
double rd= CGAL::to_double(roots[i]);
//bool ie= iem(roots[i]);
int mult=1;
//if (ie) mult=2;
for (int k=0; k<mult; ++k) {
int mm=-1;
double md=.5;
for (It c= roots_b; c != roots_e; ++c) {
if (std::abs(rd-*c) < md && taken_maple[c-roots_b]==false) {
md= std::abs(rd-*c);
mm=c-roots_b;
}
}
if (mm != -1) {
taken_maple[mm]=true;
taken_solver[i]=true;
}
}
}
bool has_error=false;
for (unsigned int i=0; i< taken_maple.size(); ++i) {
if (!taken_maple[i]) has_error=true;
}
for (unsigned int i=0; i< taken_solver.size(); ++i) {
if (!taken_solver[i]) has_error=true;
}
if (has_error) {
std::cout << std::endl;
for (unsigned int i=0; i< taken_maple.size(); ++i) {
if (!taken_maple[i])
std::cerr << "NUMERIC ISSUE Missing " << roots_b[i] << std::endl;
}
for (unsigned int i=0; i< taken_solver.size(); ++i) {
if (!taken_solver[i]) std::cerr << "NUMERIC ISSUE Extra " << roots[i] << std::endl;
}
}
std::cout << std::endl;
}
++reps;
} while (timer.time() <3 && !verbose);
timer.stop();
if (verbose) {
std::cout << "Elapsed time: " << timer.time()/reps << std::endl;
std::cout << std::endl << std::endl;
}
else {
num_roots += total_roots;
total_time += timer.time();
}
//time_estimate(q, start, CE());
}
void clear_timings() {
total_time=0;
num_roots=0;
ce_time=0;
nce=0;
}
void write_timings() {
if (!verbose) {
std::cout.setf(std::ios::fixed);
std::cout.precision(2);
std::cout << 1000000.0*total_time/static_cast<double>(num_roots) ;
if (ce_time != 0) {
std::cout << "(" << 1000000.0*ce_time/static_cast<double>(nce) << ")";
}
std::cout << " & ";
}
}
void cleaned() { {
// polynomials whose roots have been filtered based on derivitives
Fn q= cf_(-16.614835192865492,19.540659228538015,-4);
DV a; a.push_back(3.788873605);
check_polynomial(q,a.begin(), a.end(), 1.096291202);
}
{
Fn q= cf_(16.614835192865492,-19.540659228538015,4);
DV a; a.push_back(1.09630); a.push_back(3.788873605);
check_polynomial(q,a.begin(), a.end(), 1.096291202);
}
{
Fn q= cf_(3.788873605, -3.788873605);
DV a; a.push_back(1.0);
check_polynomial(q,a.begin(), a.end(), 1);
}
{
Fn q= cf_(-3.788873605, 3.788873605);
DV a;
check_polynomial(q,a.begin(), a.end(), 1);
}
{
Fn q= cf_(-2295485086.0,2072822157.0, 116461914.2, -116175036.5, -10063149.87,
-196007.0344, 3460.886000, 136.9100396, 1.0);
DV a;
a.push_back(-79.34012316);
a.push_back(-45.32403162);
a.push_back(-32.43524488);
a.push_back(-22.45230969);
a.push_back(-5.239841933);
a.push_back(1.123932833);
a.push_back(3.435646759);
a.push_back(43.32193209);
check_polynomial(q, a.begin(), a.end(), Rt(-79.34012316));
check_polynomial(q, a.begin()+2, a.end(), Rt(-45.32403162));
}
}
void wilkinson() {
clear_timings();
unsigned int max_root = 20;
CGAL_POLYNOMIAL_NS::Wilkinson_generator<K> gen(k_);
Fn q=gen(max_root);
DV a;
for (unsigned int i = 1; i <= max_root; i++) {
a.push_back(i);
}
check_polynomial(q,a.begin(), a.end());
write_timings();
}
void mignotte() {
clear_timings();
unsigned int degree = 50;
CGAL_POLYNOMIAL_NS::Mignotte_generator<K> gen(k_);
Fn q=gen(degree);
DV a;
a.push_back(-1.0925395143411487);
a.push_back(0.20000000000000001);
a.push_back(0.20000000000000001);
a.push_back(1.0756542734086822);
check_polynomial(q,a.begin(), a.end());
write_timings();
}
void square_free() {
clear_timings();
{
Fn q= cf_(1,3,-4,-7);
DV a; a.push_back(-0.8746784187); a.push_back(-0.2800178096);
a.push_back(0.5832676569);
check_polynomial(q,a.begin(), a.end());
}
{
Fn q= cf_(2,0,1);
DV a;
check_polynomial(q,a.begin(), a.end());
}
/*
{
Fn q= Fn::make(0);
DV a;
check_square_free_polynomial(q,a.begin(), a.end(),print_multiplicity);
}
{
Fn q= Fn::make(-3);
DV a;
check_square_free_polynomial(q,a.begin(), a.end(),print_multiplicity);
}
*/
{
Fn q= cf_(2,-1);
DV a;
a.push_back(2);
check_polynomial(q,a.begin(), a.end());
}
{
Fn q= cf_(1,3,-4);
DV a;a.push_back(-.25); a.push_back(1);
check_polynomial(q,a.begin(), a.end());
}
{
Fn q= cf_(2,5,-4,-7,8,-14);
DV a;a.push_back(0.7382979184);
check_polynomial(q,a.begin(), a.end());
}
{
Fn q= cf_(2, 0, -1, 1, -3, -2);
DV a;a.push_back(-1.833007477);
a.push_back(-.8474637161);
a.push_back(.7987710980);
check_polynomial(q,a.begin(), a.end());
}
write_timings();
}
void small_intervals() {
clear_timings();
Fn q= cf_(-2295485086.0,2072822157.0, 116461914.2, -116175036.5, -10063149.87,
-196007.0344, 3460.886000, 136.9100396, 1.0);
DV a;
a.push_back(-79.34012316);
a.push_back(-45.32403162);
a.push_back(-32.43524488);
a.push_back(-22.45230969);
a.push_back(-5.239841933);
a.push_back(1.123932833);
a.push_back(3.435646759);
a.push_back(43.32193209);
check_polynomial(q, a.begin()+7, a.end()-1, Rt(5.0), Rt(40.0));
check_polynomial(q, a.begin()+5, a.end(), Rt(0.0), Rt(45.0));
check_polynomial(q, a.begin()+6, a.end()-1, Rt(1.123932834), Rt(20.0));
write_timings();
}
void non_simple() {
clear_timings();
{
Fn q= cf_(-4365900, -5040, 1274073, 72979, -87647, -1045,
2275, -79, -17, 1);
// -9., -5., -3., -3., 2., 7., 7., 10., 11.
/*DV a;a.push_back(RP(-9,1));
a.push_back(RP(-5,1));
a.push_back(RP(-3,2 ));
a.push_back(RP(2,1));
a.push_back(RP(7,2));
a.push_back(RP(10,1));
a.push_back(RP(11,1));*/
DV a;a.push_back(-9);
a.push_back(-5);
a.push_back(-3);
a.push_back(-3);
a.push_back(2);
a.push_back(7);
a.push_back(7);
a.push_back(10);
a.push_back(11);
check_polynomial(q,a.begin(), a.end());
check_polynomial(q,a.begin()+2, a.end()-2, Rt(-4), Rt(9));
check_polynomial(q,a.begin()+2, a.end()-4, Rt(-4), Rt(7));
check_polynomial(q,a.begin()+2, a.end()-2, Rt(-3), Rt(9));
check_polynomial(q,a.begin()+2, a.end()-4, Rt(-3), Rt(7));
}
write_timings();
}
void roots() {
CGAL_assertion(std::numeric_limits<Rt>::has_infinity);
Rt inf= std::numeric_limits<Rt>::infinity();
Rt minf= -std::numeric_limits<Rt>::infinity();
Rt big(100000);
Rt small_rt(-1000000);
Rt z(0);
CGAL_assertion(big < inf);
CGAL_assertion(big > minf);
CGAL_assertion(small_rt < big);
CGAL_assertion(z > small_rt);
CGAL_assertion(-big < big);
CGAL_assertion(small_rt < inf);
CGAL_assertion(small_rt > minf);
CGAL_assertion(z > minf);
CGAL_assertion(z < inf);
}
void all() {
square_free();
wilkinson();
mignotte();
small_intervals();
non_simple();
}
void exact() { {
Fn q= cf_(NT("1225969745589853940699882447880155625/664613997892457936451903530140172288"),
NT("-149192919905533219658325090294590900625/2658455991569831745807614120560689152"),
NT("-176473560706138356181046066127379288725/2658455991569831745807614120560689152"),
NT("-59323957331081368745121550919931250845/10633823966279326983230456482242756608"),
NT("213976196148401335217623674372070736157/10633823966279326983230456482242756608"),
NT("10107182444809889626671506586504631263/5316911983139663491615228241121378304"),
NT("-12967429119484147645038507190192518603/10633823966279326983230456482242756608"),
NT("-894074013945673684338239220768251181/10633823966279326983230456482242756608"),
NT("963347446801307512325971073199501/5316911983139663491615228241121378304"));
DV a;
a.push_back(-14.46710600);
a.push_back(-3.440766418);
a.push_back(0.031679482287);
a.push_back(2.562471103);
a.push_back(3.490003117);
a.push_back(478.0766178);
check_polynomial(q, a.begin(), a.end());
Root_container s=
k_.root_container_object(q, Rt(-10000.0),
Rt(10000.0));
std::vector<Rt> roots;
std::copy(s.begin(), s.end(), std::back_inserter(roots));
//typename Root_container::iterator it=s.begin();
for (unsigned int offset=0; offset < roots.size()-1; ++offset) {
check_polynomial(q, a.begin()+offset, a.end(), roots[offset]);
}
}
{
Fn q= cf_(NT("2102889623460648672193194772242161145741209860309/2854495385411919762116571938898990272765493248"),
NT("-27785292706190294695336118341445556791078034998297/5708990770823839524233143877797980545530986496"),
NT("26399617884277473686036041581122386808181680932651/2854495385411919762116571938898990272765493248"),
NT("-65515073263163417483134441553521565/5070602400912917605986812821504"),
NT("6369956942652855102749740125692887/316912650057057350374175801344"),
NT("-19558908001688696225/2251799813685248"),
NT("5694238778856570877/562949953421312"),
NT("-6468/1"));
DV a;
a.push_back(0.2307182653);
a.push_back(0.4663233325);
a.push_back(1.602781661);
check_polynomial(q, a.begin(), a.end(), 0);
Root_container s=
k_.root_container_object(q, -Rt(10000.0),
Rt(100000.0));
std::vector<Rt> roots;
std::copy(s.begin(), s.end(), std::back_inserter(roots));
//typename Root_container::iterator it=s.begin();
for (unsigned int offset=0; offset < roots.size()-1; ++offset) {
check_polynomial(q, a.begin()+offset, a.end(), roots[offset]);
}
}
}
};
//typedef CGAL::Polynomial::Polynomial_kernel<Polynomial_gmpq> Kernel_gmpq;
//typedef CGAL::Polynomial::Polynomial_kernel<Polynomial_double> Kernel_double;
/*typedef CGAL::Polynomial::Filtered_kernel<Polynomial_gmpq,
Polynomial_double,
CGAL::Polynomial::Polynomial<CGAL::Interval_nt_advanced>,
Polynomial_gmpq> Filtered_kernel_gmpq;*/
//typedef CGAL::Polynomial::Polynomial_kernel<Polynomial_bigint> Kernel_core;
//typedef Kernel_gmpq Exact_kernel;
//typedef Filtered_kernel_gmpq Filtered_kernel;
//typedef Kernel_double Numeric_kernel;
#endif
Reference in New Issue View Git Blame Copy Permalink