diff options
Diffstat (limited to '')
| -rw-r--r-- | src/toys/root-finder-comparer.cpp | 247 |
1 files changed, 247 insertions, 0 deletions
diff --git a/src/toys/root-finder-comparer.cpp b/src/toys/root-finder-comparer.cpp new file mode 100644 index 0000000..3952d7c --- /dev/null +++ b/src/toys/root-finder-comparer.cpp @@ -0,0 +1,247 @@ +#include <2geom/sbasis.h> +#include <2geom/bezier-to-sbasis.h> +#include <2geom/solver.h> +#include <2geom/sbasis-poly.h> +#include "../2geom/orphan-code/nearestpoint.cpp" // FIXME: This looks like it may give problems later, (including a .cpp file) +#include <2geom/path.h> + +#include <toys/path-cairo.h> +#include <toys/toy-framework-2.h> + +#define ZROOTS_TEST 0 +#if ZROOTS_TEST +#include <2geom/zroots.c> +#endif + +using std::swap; + +namespace Geom{ +extern void subdiv_sbasis(SBasis const & s, + std::vector<double> & roots, + double left, double right); +}; + +double eval_bernstein(double* w, double t, unsigned N) { + double Vtemp[2*N]; + //const int degree = N-1; + for (unsigned i = 0; i < N; i++) + Vtemp[i] = w[i]; + + /* Triangle computation */ + const double omt = (1-t); + //Left[0] = Vtemp[0]; + //Right[degree] = Vtemp[degree]; + double *prev_row = Vtemp; + double *row = Vtemp + N; + for (unsigned i = 1; i < N; i++) { + for (unsigned j = 0; j < N - i; j++) { + row[j] = omt*prev_row[j] + t*prev_row[j+1]; + } + //Left[i] = row[0]; + //Right[degree-i] = row[degree-i]; + swap(prev_row, row); + } + return prev_row[0]; +} + +#include <vector> +using std::vector; +using namespace Geom; + +#ifdef HAVE_GSL +#include <complex> +using std::complex; +#endif + +class RootFinderComparer: public Toy { +public: + PointSetHandle psh; + + void draw(cairo_t *cr, std::ostringstream *notify, int width, int height, bool save, std::ostringstream *timer_stream) override { + std::vector<Geom::Point> trans; + trans.resize(psh.size()); + for(unsigned i = 0; i < handles.size(); i++) { + trans[i] = psh.pts[i] - Geom::Point(0, height/2); + } + + Timer tm; + + + std::vector<double> solutions; + solutions.resize(6); + + tm.ask_for_timeslice(); + tm.start(); + FindRoots(&trans[0], 5, &solutions[0], 0); + Timer::Time als_time = tm.lap(); + *notify << "original time = " << als_time << std::endl; + + D2<SBasis> test_sb = psh.asBezier();//handles_to_sbasis(handles.begin(),5); + Interval bs = *bounds_exact(test_sb[1]); + cairo_move_to(cr, test_sb[0](0), bs.min()); + cairo_line_to(cr, test_sb[0](1), bs.min()); + cairo_move_to(cr, test_sb[0](0), bs.max()); + cairo_line_to(cr, test_sb[0](1), bs.max()); + cairo_stroke(cr); + *notify << "sb bounds = "<<bs.min()<< ", " <<bs.max()<<std::endl; + Poly ply = sbasis_to_poly(test_sb[1]); + ply = Poly(height/2) - ply; + cairo_move_to(cr, 0, height/2); + cairo_line_to(cr, width, height/2); + cairo_stroke(cr); +#ifdef HAVE_GSL + vector<complex<double> > complex_solutions; + complex_solutions = solve(ply); +#if 1 + *notify << "gsl: "; + std::copy(complex_solutions.begin(), complex_solutions.end(), std::ostream_iterator<std::complex<double> >(*notify, ",\t")); + *notify << std::endl; +#endif +#endif + +#if ZROOTS_TEST + fcomplex a[ply.size()]; + for(unsigned i = 0; i < ply.size(); i++) { + a[i] = ply[i]; + } + //copy(a, a+ply.size(), ply.begin()); + fcomplex rts[ply.size()]; + + zroots(a, ply.size(), rts, true); + for(unsigned i = 0; i < ply.size(); i++) { + if(! a[i].imag()) + solutions[i] = a[i].real(); + } +#endif + +#ifdef HAVE_GSL + + tm.ask_for_timeslice(); + tm.start(); + solve(ply); + als_time = tm.lap(); + *timer_stream << "gsl poly = " << als_time << std::endl; +#endif + + #if ZROOTS_TEST + tm.ask_for_timeslice(); + tm.start(); + zroots(a, ply.size(), rts, false); + als_time = tm.lap(); + *timer_stream << "zroots poly = " << als_time << std::endl; + #endif + + #if LAGUERRE_TEST + tm.ask_for_timeslice(); + tm.start(); + Laguerre(ply); + als_time = tm.lap(); + *timer_stream << "Laguerre poly = " << als_time << std::endl; + complex_solutions = Laguerre(ply); + + #endif + + #define SBASIS_SUBDIV_TEST 0 + #if SBASIS_SUBDIV_TEST + tm.ask_for_timeslice(); + tm.start(); + subdiv_sbasis(-test_sb[1] + Linear(3*width/4), + rts, 0, 1); + als_time = tm.lap(); + *timer_stream << "sbasis subdivision = " << als_time << std::endl; + #endif + #if 0 + tm.ask_for_timeslice(); + tm.start(); + solutions.resize(0); + find_parametric_bezier_roots(&trans[0], 5, solutions, 0); + als_time = tm.lap(); + *timer_stream << "solver parametric = " << als_time << std::endl; + #endif + double ys[trans.size()]; + for(unsigned i = 0; i < trans.size(); i++) { + ys[i] = trans[i][1]; + double x = double(i)/(trans.size()-1); + x = (1-x)*height/4 + x*height*3/4; + draw_handle(cr, Geom::Point(x, height/2 + ys[i])); + } + + solutions.resize(0); + tm.ask_for_timeslice(); + tm.start(); + find_bernstein_roots(ys, 5, solutions, 0, 0, 1, false); + als_time = tm.lap(); + *notify << "found sub solutions at:\n"; + std::copy(solutions.begin(), solutions.end(), std::ostream_iterator<double >(*notify, ",")); + *notify << "solver 1d bernstein subdivision n_slns = " << solutions.size() + << ", time = " << als_time << std::endl; + + solutions.resize(0); + tm.ask_for_timeslice(); + tm.start(); + find_bernstein_roots(ys, 5, solutions, 0, 0, 1, true); + als_time = tm.lap(); + + *notify << "solver 1d bernstein secant subdivision slns" << solutions.size() + << ", time = " << als_time << std::endl; + *notify << "found secant solutions at:\n"; + std::copy(solutions.begin(), solutions.end(), std::ostream_iterator<double >(*notify, ",")); + *notify << "solver 1d bernstein subdivision accuracy:" + << std::endl; + for(double solution : solutions) { + *notify << solution << ":" << eval_bernstein(ys, solution, trans.size()) << ","; + } + tm.ask_for_timeslice(); + tm.start(); + solutions = roots( -test_sb[1] + Linear(height/2)); + als_time = tm.lap(); +#if 1 + *notify << "sbasis roots: "; + std::copy(solutions.begin(), solutions.end(), std::ostream_iterator<double>(*notify, ",\t")); + *notify << "\n time = " << als_time << std::endl; +#endif + for(double solution : solutions) { + double x = test_sb[0](solution); + draw_cross(cr, Geom::Point(x, height/2)); + + } + + *notify << "found " << solutions.size() << "solutions at:\n"; + std::copy(solutions.begin(), solutions.end(), std::ostream_iterator<double >(*notify, ",")); + + D2<SBasis> B = psh.asBezier();//handles_to_sbasis(handles.begin(), 5); + Geom::Path pb; + pb.append(B); + pb.close(false); + cairo_path(cr, pb); + + B[0] = Linear(width/4, 3*width/4); + cairo_d2_sb(cr, B); + Toy::draw(cr, notify, width, height, save,timer_stream); + } + RootFinderComparer(unsigned degree) + { + for(unsigned i = 0; i < degree; i++) psh.push_back(Geom::Point(uniform()*400, uniform()*400)); + handles.push_back(&psh); + } +}; + +int main(int argc, char **argv) { + unsigned bez_ord = 6; + if(argc > 1) + sscanf(argv[1], "%d", &bez_ord); + init(argc, argv, new RootFinderComparer(bez_ord)); + + return 0; +} + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 : |