
12

I am stuck. I got an implementation of quicksort as an example of list comprehension in the manual:
function quicksort(arr) =
!(len(arr)>0) ? [] :
let( pivot = arr[floor(len(arr)/2)],
lesser = [ for (y = arr) if (y < pivot) y ],
equal = [ for (y = arr) if (y == pivot) y ],
greater = [ for (y = arr) if (y > pivot) y ]
)
concat( quicksort(lesser), equal, quicksort(greater) );
It works well for very large random arrays and, in the worst case, for arrays with length at least 3000 without stack overflow.
Usually the partitions of quicksort are in two subarrays. So I tweaked the code to get:
function quicksort1(arr) =
!(len(arr)>0) ? [] :
let( pivot = arr[floor(len(arr)/2)],
lesser = [ for (y = arr) if (y <= pivot) y ],
greater = [ for (y = arr) if (y > pivot) y ]
)
concat( quicksort1(lesser), quicksort1(greater) );
where the partition "equal" is included in the partition "lesser".
The trouble is that this last version is not working: I get a stack overflow for any non void array. And I don't see why.


In the first implementation there is always at least one element that
gets taken away from the sorting (the pivot), so that you will always
finish.
In quicksort1 when the array has just 1 element it gets always selected
as pivot, always included in lesser and you get infinite recursion.
> I am stuck. I got an implementation of quicksort as an example of
> list
> comprehension in the manual:
>
> > function quicksort(arr) =
> > !(len(arr)>0) ? [] :
> > let( pivot = arr[floor(len(arr)/2)],
> > lesser = [ for (y = arr) if (y < pivot) y ],
> > equal = [ for (y = arr) if (y == pivot) y ],
> > greater = [ for (y = arr) if (y > pivot) y ]
> > )
> > concat( quicksort(lesser), equal, quicksort(greater) );
>
> It works well for very large random arrays and, in the worst case, for
> arrays with length at least 3000 without stack overflow.
> Usually the partitions of quicksort are in two subarrays. So I
> tweaked the code to get:
>
> > function quicksort1(arr) =
> > !(len(arr)>0) ? [] :
> > let( pivot = arr[floor(len(arr)/2)],
> > lesser = [ for (y = arr) if (y <= pivot) y ],
> > greater = [ for (y = arr) if (y > pivot) y ]
> > )
> > concat( quicksort1(lesser), quicksort1(greater) );
>
> where the partition "equal" is included in the partition "lesser".
> The trouble is that this last version is not working: I get a stack
> overflow for any non void array. And I don't see why.
>
>
>
> 
> View this message in context:
> http://forum.openscad.org/Quicksorttp17044.html Sent from the
> OpenSCAD mailing list archive at Nabble.com.
>
> _______________________________________________
> OpenSCAD mailing list
> [hidden email]
> http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


Bingo! Why I didn't see that? Well, I think I hadn't fully grasped quicksort algorithm. :(
Meanwhile, I got a version of quicksort which seems to have the pattern of tail recursion:
function qs(arr, arr_out=[]) =
!(len(arr)>0) ?
arr_out :
qs( arr = high_part(arr),
arr_out = concat(qs( low_part(arr), arr_out ), middle_part(arr)) );
function high_part (arr) = [ for (y = arr) if (y > pivot(arr)) y ];
function middle_part(arr) = [ for (y = arr) if (y == pivot(arr)) y ];
function low_part (arr) = [ for (y = arr) if (y < pivot(arr)) y ];
function pivot (arr) = arr[floor(len(arr)/2)];
I were able to sort a worst case array of length 5000 in 48sec with qs(). The previous quicksort() function fail with stack overflow for lengths lesser than 4000. However, it is much slower: for a random array with 40000 values, qs() required 60 sec while quicksort() spent just 4sec.


Interesting.
With 40000 random numbers I also get 4s with quicksort but qs takes 3 minutes 37 seconds for me. I translated it into Python to try to see what the difference is, OpenScad desperately needs echo() in functions! import random import time
def pivot(arr): return arr[int(len(arr)/2)]
def high_part(arr): p = pivot(arr) return [y for y in arr if y > p]
def middle_part(arr): p = pivot(arr) return [y for y in arr if y == p]
def low_part(arr): p = pivot(arr) return [y for y in arr if y < p]
def qs(arr, arr_out=[], n = 0): #print n, arr return arr_out if len(arr) == 0 else qs(high_part(arr), qs(low_part(arr), arr_out, 2) + middle_part(arr), 1)
def quicksort(arr, n = 0): #print n, arr if len(arr) == 0: return [] pivot = arr[int(len(arr)/2)] lesser = [y for y in arr if y < pivot] equal = [y for y in arr if y == pivot] greater = [y for y in arr if y > pivot] return quicksort(lesser, 2) + equal + quicksort(greater, 1)
list = [random.random() for x in xrange(40000)]
t0 = time.time() qs(list) t1 = time.time() quicksort(list) t2 = time.time()
print "qs:", round(t1  t0, 3), "quicksort:", round(t2  t1, 3)
There is also a big difference in timings: qs: 5.632 quicksort: 0.256 I hoisted the pivot calculations out of the loops and I did the same in the openscad version with let but it made no difference.
If I add the commented out print statements and use a small array of integers it prints exactly the same calls and intermediate results, so I can't explain the big time difference. Perhaps there are more list copies due to passing arr_out forward.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


In nophead's benchmark, Python is 15.6 times faster for quicksort, and 38.5 times faster for qs.
I attribute this to the fact that Python compiles to bytecode. We should be able to get a 10x speedup pretty easily with a bytecode compiler, even before doing much performance tuning. Note that Python is slow compared to Javascript, which has insane amounts of performance optimization.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


I dunno, why would anyone want to triangulate a polygon? I've been thinking about performance issues while reading that thread. I agree that triangulating a polygon is not the sort of task that OpenSCAD was designed for.
I think that some of this is driven by the fact that OpenSCAD provides a limited set of primitive shapes, and if you want to go beyond that, you need to use polyhedron(). And if you really push the limits of what can be expressed with polyhedron(), then you can get into requiring these complex calculations.
I'd like to provide an easier way than polyhedron() for constructing complex mathematically defined shapes. I think Functional Representation would be a great thing to support. Of course, that opens up an even bigger can of worms than simply speeding up the interpreter using a byte code compiler.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


Only Ronaldo can explain why one would need quicksort and polygon triangulation in OpenSCAD. He started both of those threads. Maybe he will explain.
You might need 40,000 vertices in a polygon as an end cap for a tube, either a really accurate circle, or some other shape. Fractal geometry is another general reason for needing lots of vertices.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


As far as I know, polyhedron faces isn't triangulated in rendering just in preview. That is the main reason I searched for a polygon triangulation. Apparently, when CGAL rejects a face as nonplanar, OpenSCAD tries to triangulate it: the alternative procedure. But sometimes this doesn't work either.
I don't need to sort 40000 array values. I was only doing a time analysis to compare different methods. Remember that the quicksort exemplified in the manual, to be robust, should be restricted to handle vectors with 3000 elements. That is the limit for the worst (very unlikely) case, but you never know.


As I said before in another thread, OpenSCAD should provide a good set of list processing operators. And sorting is a very basic one. After all, lists are the only available data structure in the language.
Besides, CGAL has lots of computational geometry methods and very very few of them are accessible by the language. My dream is to have a general OpenSCAD interface to CGAL library for people who wants to dig in the CGAL world.
To dough: I like Frep too. It is a very general paradigm but it requires some mathematicaloriented mind. I see no way to include Frep without a robust optimized kernel to support it. If CGAL is able to give that support, it would a good path to pursuit.


As I wrote in the other thread, quicksort alone will not speed up your code from O(n²) to O(n log n). What you need, is a sorter operating over a balanced tree. Thus: better "ask" for this. BUT: this approach needs function pointers, as you will have to provide a predicate for general sorting. So better wait for OpenSCAD2, use a richer programming language, or try your luck and ask for an API ;)


Ronaldo wrote: As far as I know, polyhedron faces isn't triangulated in rendering just in preview. That is the main reason I searched for a polygon triangulation. Apparently, when CGAL rejects a face as nonplanar, OpenSCAD tries to triangulate it: the alternative procedure. But sometimes this doesn't work either.
Ronaldo, it sounds like you are writing your own polygon triangulation software in order to work around a possible bug in polyhedron(). There is some argument to polyhedron(), involving a nonplanar face, that works in preview, but doesn't work in rendering.
Can you describe the bug in more detail, maybe provide a failing test case?
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


doug.moen wrote
Ronaldo, it sounds like you are writing your own polygon triangulation
software in order to work around a possible bug in polyhedron(). There is
some argument to polyhedron(), involving a nonplanar face, that works in
preview, but doesn't work in rendering.
Can you describe the bug in more detail, maybe provide a failing test case?
I don't have a failing test case anymore. I integrated the polygon triangulation in my system and have no more the issues I had before.
I am working on a system that allows solid (sets with manifold boundaries) modelling using polyhedron primitive with faces defined by Bezier patches and planar (or almost) faces. It evolved very well but I faced now and then with annoying messages by CGAL of nonplanar faces. The analysis of the face definition usually shown that the face was planar (for instance, all vertices had the same x coordinate). Sometimes the alternate construction solved the trouble. Sometimes not. So, I decided to write my own polygon triangulation algorithm which unhappily is not good enough and took too much of my time. The discussions on simple polygon test, plane fitting and sorting, all came from that target. That is the story.
I have no way to say that I was facing a bug. I don't know what is being considered for OpenSCAD2. But, it would be nice to have access to this sort of geometrical methods that are in the core of CGAL. And a richer set of list processing functions.

Administrator

> On Apr 12, 2016, at 11:01 AM, Ronaldo < [hidden email]> wrote:
>
> Besides, CGAL has lots of computational geometry methods and very very few
> of them are accessible by the language. My dream is to have a general
> OpenSCAD interface to CGAL library for people who wants to dig in the CGAL
> world.
>
Ronaldo,
Keep in mind that OpenSCAD isn’t primarily a wrapper around CGAL. We use it as a last resort to perform certain types of computational geometry operations. Due to severe performance challenges with CGAL, we’re actively looking for ways to reduce our dependency on CGAL.
To dig into the CGAL world, I would recommend using their officially supported SWIG bindings as a starting point.
Marius
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org

Administrator

> On Apr 12, 2016, at 10:44 AM, Ronaldo < [hidden email]> wrote:
>
> Apparently, when CGAL rejects a face as nonplanar, OpenSCAD tries to
> triangulate it: the alternative procedure. But sometimes this doesn't work
> either.
>
“sometimes this doesn’t work” sounds like a bug.
CGAL frequently rejects polygons as nonplanar, if those polygons were specified using floating point vertices. We have workarounds for that, but our workaround are only as good as our test cases. Any additional corner cases breaking this would be valuable input allowing us to be more robust.
As a sidenote for anyone wondering: The polygon triangulator built into CGAL isn’t very robust and crashes on some of our test cases. We thus had to switch to a more robust implementation, currently libtess2, as mentioned earlier in this thread.
Marius
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


"CGAL frequently rejects polygons as nonplanar, if those polygons were specified using floating point vertices. We have workarounds for that, but our workaround are only as good as our test cases."
I suggest that polyhedron() should automatically triangulate any face with more than 3 vertices, since in this case floating point vertices are always used.
_______________________________________________
OpenSCAD mailing list
[hidden email]
http://lists.openscad.org/mailman/listinfo/discuss_lists.openscad.org


kintel wrote
Keep in mind that OpenSCAD isn’t primarily a wrapper around CGAL. We use it as a last resort to perform certain types of computational geometry operations. Due to severe performance challenges with CGAL, we’re actively looking for ways to reduce our dependency on CGAL.
To dig into the CGAL world, I would recommend using their officially supported SWIG bindings as a starting point.
Thank you for the clarification and reference. I imagined that. And it was the reason why I did not submitted a feature request following a Parkinbot suggestion. Anyway it was just a dream...:)

12
