Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

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 :