Hello,
since some days I ask myself if it is already possible with Nimrod to write an elegant convex hull algorithm similar to this wikipedia page:
http://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain
(I wrote the Ruby code some years ago, based on the Python example...)
The algorithm itself should be no problem, but generally we remove duplicated points at the beginning -- is there something like uniq! already available for nimrods's arrays. Or have we to push all array elements to a hash manually, and back again? I was not able to find somethings like an overview of array methods for nimrod at all, only array sorting.
And finally, the Python and Ruby examples have a one or two line minimal test at the end. Currently I have not really an idea for such a compact test in Nimrod.
Best regards,
Stefan Salewski
You can use sequtils.distnct (misspelling intentional) or toSeq(toSet(a).items) with sequtils.toSeq and sets.toSet to remove duplicates, but it's faster to make use of the fact that the list is already sorted, if your goal isn't code golf (short over fast code), especially since sequtils.distnct currently has quadratic time complexity:
import algorithm
type Point = tuple[x, y: float]
proc cmpPoint(a, b: Point): int =
result = cmp(a.x, b.y)
if result == 0:
result = cmp(a.y, b.y)
proc sortUniq[T](a: var seq[T], cmp: proc(x, y: T): int {.closure.}) =
if len(a) > 0:
a.sort(cmp)
var p = 0
for i in 1..high(a):
if a[i] != a[p]:
p += 1
a[p] = a[i]
a.setLen(p+1)
You can generate the desired list for a quick test as follows:
import sequtils
echo map(toSeq(0..99), proc(x: int): Point = (float(x div 10), float(x mod 10)))
Or, using the new lambda syntax:
import sequtils, future
echo map(toSeq(0..99), (x: int) => (float(x div 10), float(x mod 10)))
seqs can be compared for equality using ==.
import algorithm, sequtils
type Point = tuple[x, y: float]
proc cmpPoint(a, b: Point): int =
result = cmp(a.x, b.x)
if result == 0:
result = cmp(a.y, b.y)
proc sortUniq[T](a: var seq[T], cmp: proc(x, y: T): int {.closure.}) =
#if len(a) > 0:
a.sort(cmp)
var p = 0
for i in 1..high(a):
if a[i] != a[p]:
p += 1
a[p] = a[i]
a.setLen(p + 1)
template cross[T](o, a, b: T): expr =
(a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x)
proc convex_hull[T](points: var seq[T], cmp: proc(x, y: T): int {.closure.}) : seq[T] =
if len(points) < 2: return points
sortUniq(points, cmp)
if len(points) < 3: return points
var lower: seq[T] = @[]
for p in points:
while len(lower) > 1 and cross(lower[len(lower) - 2], lower[len(lower) - 1], p) <= 0:
discard lower.pop
lower.add(p)
var upper: seq[T] = @[]
for i in countdown(high(points), low(points)):
let p = points[i]
while len(upper) > 1 and cross(upper[len(upper) - 2], upper[len(upper) - 1], p) <= 0:
discard upper.pop
upper.add(p)
result = concat(lower[0 .. -2], upper[0 .. -2])
var s = map(toSeq(0..99), proc(x: int): Point = (float(x div 10), float(x mod 10)))
assert(convex_hull[Point](s, cmpPoint) == @[(0.0, 0.0), (9.0, 0.0), (9.0, 9.0), (0.0, 9.0)])Negative indices to index an array from the end cannot work for general arrays because array indices can be negative normally:
var s: array[-1..1, int]
The usual idiom to index something from the end is to use a[a.high-k].
With respect to initializing sequence variables, you can also use
var lower = newSeq[T]()
to initialize a seq variable with an empty list of a specific type.
Negative indices to index an array from the end cannot work
OK, I missed that. It is interesting that we now have arbitrary ranges for arrays again -- I think we had it already in Pascal, but Wirth forces arrays to start at index 0 for Oberon, as in C. His motivation was, that start at 0 is fine in most cases, that other start index can lead to confusion, and that offset slow down and increase code size. I really loved the negative indices in Ruby...
var lower = newSeqT
Fine. I read about newSeq in tutorial or manual, but did not know how to use it.
I'm not sure that saving the effort to type a.high is worth that.
I think accessing an array in reverse order is not a rare task. Generally a backward iterator will do.
Accessing the last element with an expression like hull[hull.high] compared to hull[-1] is not really nice -- for writing and reading. Maybe something like hull.last may be possible?
I discovered the negative indices some years ago when learning Ruby -- first using it was a bit strange, but it makes sense when we regard array elements as a circular list -- predecessor of index 0 is -1. Performance is indeed a good reason against it.
Stefan: I think accessing an array in reverse order is not a rare task. Generally a backward iterator will do.
Negative indices don't help you here, since you can't do countdown(-1, 0) and have it work (because -1 < 0). You'll have to do something like:
for i in countdown(a.high, a.low):
doSomethingWith(a[i])
Stefan: Accessing the last element with an expression like hull[hull.high] compared to hull[-1] is not really nice -- for writing and reading. Maybe something like hull.last may be possible?
An a.last function in the standard library (possibly along with a.first) would indeed be nice. I have one in my own odds-and-ends library, along with a bunch of other helpful array/seq functionality, but it's something that should probably be on a "batteries included" list.
Stefan: I discovered the negative indices some years ago when learning Ruby -- first using it was a bit strange, but it makes sense when we regard array elements as a circular list -- predecessor of index 0 is -1. Performance is indeed a good reason against it.
Performance, oddly enough, is one of the major reasons why you commonly find such a construct in interpreted languages. In compiled languages, doing -1 over len(a)-1 does not buy you anything, since the runtime has to insert a separate branch and the additional computation, but in an interpreted language, having the interpreter do the extra len(a) and subtraction can be quite expensive for an inner loop.
This is also why it's generally not as big a deal if such functionality doesn't exist in the standard library: Unlike with interpreters, you do not have to fiddle with the internals of the system, but can simply write your own function (in Nimrod, a proc or template) to implement it. E.g.:
proc last*[T](a: openarray[T]): T {.inline.} = a[a.high]
iterator reverseItems[T](a: openarray[T]): T =
for i in countdown(a.high, 0):
yield a[i]
Below is my first attempt of a parallel computation of the convex full, parallel for lower and upper part. Indeed it works very well, even with generic types. Unfortunately the parallel code below takes nearly the same amount of time as the plain code above -- on my two core AMD64 box. But I really think that is not bad, the algorithm is not really good for parallel computation. The Python code from the wikibook page is about 100 times slower with Python2 (and fails with Python3).
EDIT:
I have just discovered that generation of the initial sequence of points (var s = map() below) takes already a large part of the total time. But still my feeling is that the parallel algorithm is a bit slower, tested with unix time command.
EDIT2:
I have to admit that the Nim code is only 40 times faster than the Python2 from wikibook. Reason: For my final tests I compiled with -d:release and ignored that this removes the assert statements, so an echo operation is necessary to ensure that all code is really executed.
# parallel convex hull for Nim bigbreak
# nim c --threads:on -d:release pconvex_hull.nim
import algorithm, sequtils, threadpool
type Point = tuple[x, y: float]
proc cmpPoint(a, b: Point): int =
result = cmp(a.x, b.x)
if result == 0:
result = cmp(a.y, b.y)
proc sortUniq[T](a: var seq[T], cmp: proc(x, y: T): int {.closure.}) =
a.sort(cmp)
var p = 0
for i in 1..high(a):
if a[i] != a[p]:
p += 1
a[p] = a[i]
a.setLen(p + 1)
template cross[T](o, a, b: T): expr =
(a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x)
proc half[T](p: seq[T]; upper: bool): seq[T] =
var h: seq[T] = @[]
var d, i: int
if upper:
d = 1
i = 0
else:
d = -1
i = high(p)
while i >= low(p) and i <= high(p):
while len(h) > 1 and cross(h[len(h) - 2], h[len(h) - 1], p[i]) <= 0:
discard h.pop
h.add(p[i])
i += d
discard h.pop
return h
proc convex_hull[T](points: var seq[T], cmp: proc(x, y: T): int {.closure.}) : seq[T] =
if len(points) < 2: return points
sortUniq(points, cmp)
if len(points) < 3: return points
var ul: array[2, FlowVar[seq[T]]]
parallel:
for k in 0..ul.high:
ul[k] = spawn half[T](points, k == 0)
result = concat(^ul[0], ^ul[1])
var s = map(toSeq(0..999999), proc(x: int): Point = (float(x div 1000), float(x mod 1000)))
assert(convex_hull[Point](s, cmpPoint) == @[(0.0, 0.0), (999.0, 0.0), (999.0, 999.0), (0.0, 999.0)])I would just skip duplicate points
Indeed -- already did it during an attempt to optimize the code. The result is below, still a bit slower than the initial code without parallel processing.
# parallel convex hull for Nim bigbreak
# nim c --threads:on -d:release pconvex_hull.nim
import algorithm, sequtils, threadpool
type Point = tuple[x, y: float]
proc cmpPoint(a, b: Point): int =
result = cmp(a.x, b.x)
if result == 0:
result = cmp(a.y, b.y)
template cross[T](o, a, b: T): expr =
(a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x)
template pro(): expr =
while lr1 > 0 and cross(result[lr1 - 1], result[lr1], p[i]) <= 0:
discard result.pop
lr1 -= 1
result.add(p[i])
lr1 += 1
proc half[T](p: seq[T]; upper: bool): seq[T] =
var i, lr1: int
result = @[]
lr1 = -1
if upper:
i = 0
while i <= high(p):
pro()
i += 1
else:
i = high(p)
while i >= low(p):
pro()
i -= 1
discard result.pop
proc convex_hull[T](points: var seq[T], cmp: proc(x, y: T): int {.closure.}) : seq[T] =
if len(points) < 2: return points
points.sort(cmp)
var ul: array[2, FlowVar[seq[T]]]
parallel:
for k in 0..ul.high:
ul[k] = spawn half[T](points, k == 0)
result = concat(^ul[0], ^ul[1])
var s = map(toSeq(0..999999), proc(x: int): Point = (float(x div 1000), float(x mod 1000)))
echo convex_hull[Point](s, cmpPoint)
assert(convex_hull[Point](s, cmpPoint) == @[(0.0, 0.0), (999.0, 0.0), (999.0, 999.0), (0.0, 999.0)])The main reason why you don't get a speedup is the program spends most of its time copying a sequence of 1,000,000 tuples of floats from one thread to another. Message passing is generally not a good idea when you're dealing with such large amounts of data.
Here is a screenshot of profiling results for your code (run 10x). You'll see that the half procedure accounts for only a small fraction of overall execution time. If you want to speed this up for large data, the data needs to be shared between threads (so as to avoid copying).
just added a 3d convex hull pure nim implementation w/ adjacent faces, perfect for waterman polyhedron drawing, includes a opengl interface
But where is the link to your own published research paper? Or links to used papers?
It is funny that you used the famous halve edge data structure. I used that one too for my fully dynamic constrained delaunay triangulation, but initially that halve edge ops damaged my brain.