Seems like a very good ideal overall. I get substantial speedup for the sample script with overlap=false;
3.5s vs 82s with fast-union off, about 23x speedup.

And just for fun, I wrote a variation of your script which is more cache-friendly (it makes the inner loop / entire columns of spheres cache-able)

union() {
  for (i=[0:N-1]) translate([i,0,0])
    for(j=[0:N-1]) translate([0,j,0])
      sphere(d=(overlap ? 1.1 : 0.9), $fn=50);
}

That gets my render time down to 0.732s, which is outstanding!

...However, going back to your original script and changing to overlap=true;, my testing actually shows that fast-union is significantly slower than before:
fast-union OFF: Total rendering time: 0:01:30.999
fast-union ON: Total rendering time: 0:03:08.003
So just over 2x slower with fast-union enabled.
And I similar speeds of ~1m30s on latest nightly snapshot, and even 2019.05 release as well.

// Rendering is still faster than before thanks to the caching of the spheres,

Since I can't reproduce any speedup, I'm curious what you mean about caching spheres? I haven't had a chance to look over the whole PR yet, but spheres should have already gone into Geometry Cache, no different from other Polysets.

And on the topic of prioritizing order of operations, I don't know if you've seen it, but the discussion in #1234 and my subsequent PR #3121 are maybe worth looking at.

As a quick example of why combining least-complex-geometries first, is a reasonable heuristic:
Say we have 8 objects A,B,C...H, which are all disjoint and have the same complexity.
The main cost here is basically memory accesses to read input operands, and write the output.
Since these all have same complexity, lets call it 1 Standardized Memory Access Cost Unit (SMACU)

If we just union them in order: ((((((A+B)+C)+D)+E)+F)+G)+H
Then each subsequent operation costs increasingly more:
1+1 (read) = 2 (write) (4 total SMACUs for this op)
2+1 = 3 (6 SMACUs)
3+1 = 4 (8 SMACUs)
...
7+1 = 8 (16 SMACUs)
In the end we get (4+6+8+10+12+14+16) = Grand Total 70 SMACUs

If we instead pick the lowest complexity pair between inputs and intermediate results at each step, then we end up looking more like this: ((A+B)+(C+D))+((E+F)+(G+H))
1+1=2 (4)
1+1=2 (4)
1+1=2 (4)
1+1=2 (4)
2+2 = 4 (8)
2+2 = 4 (8)
4+4 = 8 (16)
(4+4+4+4+8+8+16) = Grand Total 48 SMACUs

Differences would be even more dramatic if for example the first object A had a cost of 100, while the remaining objects stayed at cost = 1.

Anyways, using the "Morton" Z-Order curve you mentioned (i had to look it up) seems promising as well. I had also thought about mapping bounds or centroid to a position on a 3D Hilbert Curve to facilitate sorting by proximity. More or less the same idea, just using a different space filling curve.

Seems like a very good ideal overall. I get substantial speedup [...]

Thanks for trying things out!

variation of your script which is more cache-friendly

Hah, neat trick!

what you mean about caching spheres? [...] spheres should have already gone into Geometry Cache, no different from other Polysets.

Sorry I was a bit cryptic: the spheres' PolySets are indeed in the geometry cache, but when unioning N of them we convert each one to a Nef polyhedron independently (createNefPolyhedronFromGeometry in applyUnion3D) and don't cache that Nef. When I ran with $fn=100 (that is where I saw the speedup), each conversion alone was taking 7 seconds on my laptop.

The questionably-named LazyGeometry helper provides a transparent on-demand conversion & cache behaviour and I think could potentially be used beyond the union case.

I don't know if you've seen it, but the discussion in #1234 and my subsequent PR #3121 are maybe worth looking at.

I hadn't, thanks for the great pointers! I hadn't understood the extent of the problem that queue was solving, thanks for elaborating.

I'm working on refactoring the code to preserve the queue logic ("almost" done), will ping this thread when it's ready (moving the BoundingBoxes object inside the LazyGeometry, lots of the code should stay the same).

pick the lowest complexity pair between inputs and intermediate results at each step

I'd assume the most efficient way to multithread this would be work-stealing from the same queue with a proper lock?

mapping bounds or centroid to a position on a 3D Hilbert Curve to facilitate sorting by proximity.

I haven't played with Hilbert's curve yet, could be a better pick (thanks for adding the links ;-)). In any case, overlaps of the 1D intervals corresponding to a bounding box in the curve space may not result in actual overlaps of the 3D bounding boxes, so will still need some tests on the bboxes but it should make it more likely for the heuristic to find good non-overlaps in some bounded time. Can't wait to give it a try :-)

more cache-friendly (it makes the inner loop / entire columns of spheres cache-able)

Oh, and this opens an entirely new can of worms. A rendering tree optimizer might be able to detect repeating subpatterns and reproduce this manual optimization... Something I'd love to explore (see my ramblings in #3637 about tree rewrites).

I've retrofitted the original queue logic and did see sensible improvements, although haven't done thorough benchmarks yet.

(Also, I tried caching intermediate results to see if it could somehow approximate your manual optimization: it resulted in lots of cache evictions, so not a free pass just yet; left it out of the last batch of commits)

This pull request introduces 1 alert when merging 8ea55bc into 99cae25 - view on LGTM.com

new alerts:

• 1 for Inconsistent definition of copy constructor and assignment ('Rule of Two')

what you mean about caching spheres? [...] spheres should have already gone into Geometry Cache, no different from other Polysets.

Sorry I was a bit cryptic: the spheres' PolySets are indeed in the geometry cache, but when unioning N of them we convert each one to a Nef polyhedron independently (createNefPolyhedronFromGeometry in applyUnion3D) and don't cache that Nef. When I ran with $fn=100 (that is where I saw the speedup), each conversion alone was taking 7 seconds on my laptop.

I tried again with $fn=100 and i still get 2x slowdown:

N=5, overlap=true, $fn=100
nightly:    Total rendering time: 0:06:31.726
fast-union: Total rendering time: 0:12:57.206

Edit: Took me a while from starting this reply to finishing it, and in the meantime I pulled your latest changes and retested.
I am seeing some speedup now even with overlap = true;, so that's good.

N=5, overlap=true, $fn=50
nightly:    Total rendering time: 0:01:32.899
fast-union  Total rendering time: 0:01:17.869

The questionably-named LazyGeometry helper provides a transparent on-demand conversion & cache behaviour and I think could potentially be used beyond the union case.

OK, this sounds a bit similar to the smartCache* functions in GeometryEvaluator. It is supposed to handle getting from and inserting to the proper cache: GeometryCache vs CGALCache, but I think in many ways it is suboptimal, and not particularly "smart" in how it operates. Maybe the LazyGeometry functions can be used to replace this instead?
I definitely would prefer more general solutions that can be applied in multiple places over having tons of niche classes for every little thing.

I wonder about the actual complexity of Nef unions, I feel like there might be some log complexity lookups (do they use some bounding box hierarchy of faces to compute intersections?) that would make the costs even less linear?

I don't know precisely, but I'm pretty sure its very close to linear time by total number of faces. IIRC @rcolyer helped show this to me a while ago by writing some scripts and benching them at various complexities, maybe around when I was writing PR 3121.

pick the lowest complexity pair between inputs and intermediate results at each step

I'd assume the most efficient way to multithread this would be work-stealing from the same queue with a proper lock?

I don't know much about work stealing, but i guess that implies each thread would have maintain its own separate queue of work to be stolen from? I think it would be better to have a centralized queue so that ordering can be applied across threads. So I would probably start with a basic design where all threads share a single priority queue and just pull from and push to that (with locking), something like:

1. lock priority queue
2. if worker has a result (from previous loop iteration) then insert onto the queue
3. if (queue.size < 2) { release lock; break; }
4. pop next two objects
5. release lock
6. do work
7. loop to 1)

But yes proper locking will be crucial for any shared resources (and there's a lot of static / global data around the codebase, for various geometry caches, stat cache, module cache, stackchecker, progress bar, console output handlers etc. which will require careful consideration). There has been some talk, about working towards minimizing such global state wherever possible so that its more compatible with multi-threading.

I haven't played with Hilbert's curve yet, could be a better pick (thanks for adding the links ;-)). In any case, overlaps of the 1D intervals corresponding to a bounding box in the curve space may not result in actual overlaps of the 3D bounding boxes, so will still need some tests on the bboxes but it should make it more likely for the heuristic to find good non-overlaps in some bounded time. Can't wait to give it a try :-)

Another idea is maybe to just calculate volume of bbox intersections.
If you assume uniformly distributed vertices/faces, then you could estimate that the number of faces removed by a union would be roughly:
NumFaces(A) * Volume( bbox(A) ⋂ bbox(B) ) / Volume( bbox(A) ) + NumFaces(B) * Volume( bbox(A) ⋂ bbox(B) ) / Volume( bbox(B) )

I don't think this would be compatible in any way with the priority queue as-is, since that only calculates cost of individual geometries and not pairs of them.
Evaluating all possible pairs would also be a O(n^2) operation for n geometries, which might be too much overhead to do in-between every operation, maybe there could be some threshold of face count, below which the code doesn't bother to find optimal bounding box pairs etc, and falls back to priority queue. For example if a script is only doing union on a huge number axis aligned cubes, I would think checking bounding boxes etc could be as much overhead or more compared to just doing the operation itself.

Its hard to predict the impact (positive or negative) of various optimization attempts without just implementing and measuring them directly.

Oh, and this opens an entirely new can of worms. A rendering tree optimizer might be able to detect repeating subpatterns and reproduce this manual optimization...
(Also, I tried caching intermediate results to see if it could somehow approximate your manual optimization: it resulted in lots of cache evictions, so not a free pass just yet; left it out of the last batch of commits)

I don't think this type of optimization will be possible, due to the way caching works. Cache keys are basically just the substring / subtree of the CSG Tree which you see from "Design -> Display CSG Tree..."

The manual optimization works by splitting the translate into two separate operations. So the inner loop gets repeated verbatim, as a child of the outer translate. Without splitting up the translates, then each one will always be unique, so it would be pretty tough I think to match patterns like that, and even if something was devised it would probably have to be specific to grid-like translations, which is a niche solution.

Also, I will probably expand some more in your other PR eventually, but keep in mind that CSGTreeNormalizer already basically converts to a sum of products form. You can see CSGTreeNormalizer::match_and_replace for the replacement rules etc.
But this is only used for preview and not tied to CGAL evaluation. The OpenCSG rendering library requires that the CSG terms be in this form, but I have doubts that it would be helpful for CGAL due to the way that the number of terms tend to multiply to quite large values for certain combinations of ops.

I am seeing some speedup now

Yeah for me it only seems significant with large $fn values (i.e. when nefs are costly to create)

The questionably-named LazyGeometry helper provides a transparent on-demand conversion & cache behaviour and I think could potentially be used beyond the union case.

OK, this sounds a bit similar to the smartCache* functions in GeometryEvaluator. It is supposed to handle getting from and inserting to the proper cache: GeometryCache vs CGALCache, but I think in many ways it is suboptimal, and not particularly "smart" in how it operates. Maybe the LazyGeometry functions can be used to replace this instead? I definitely would prefer more general solutions that can be applied in multiple places over having tons of niche classes for every little thing.

I've seen it, and agree! I've split this logic to a mixed_cache.h file as a starting point (with a meta-templated getOrSet that picks the right cache(s) for the asked type). Would it be fine if I attempted to unify cache access in a subsequent PR? (this one is already getting large)

So I would probably start with a basic design where all threads share a single priority queue and just pull from and push to that (with locking), something like [...]

Glad we're on the same page on this.

I've actually taken a stab at this a few weeks ago: parallel_reduce.h. I wasn't using the priority queue mechanism yet, and the rest of that branch would need rethink or refresh. It does naive parallelism using futures: I essentially swapped shared_ptr<const Geometry> for lazy_ptr<const Geometry>, which is a drop-in replacement that blocks when the actual geom is dereferenced (lazy_ptr.h, pardon my clunky C++). I unlocked a few blocking calls (i.e. in caching code or everywhere that wants to know about the actual geometry's class or dimension) and make a few calls async (transforms, operations... which use that parallel_reduce, but really should probably just use TBB), and despite being globally thread-unsafe, it still seemed very stable (prolly didn't try to crash it hard enough).

My take away from that experiment was that unions were still too slow even with multithreading (and the evaluation tree is too deep, which serializes lots of operations), hence the lazy union extravaganza from #3637 (to flatten the tree) and now this fast union route.

FYI here's the real-world model that motivated this month of hacking: it renders 20x faster with all my (monothreaded) features combined.

BTW: would you be open to having some benchmark model collection + scripts? (I'm sure there's cool tools to use?) I'd be happy to contribute.

[...] calculate volume of bbox intersections. [...]

Interesting idea! (even if we don't apply it here, could help with the cost function of some global optimizer... thinking out loud, but imagine if we could build a gross, octree-based approximation of the model with just the leaf #facet & bbox information? we might have wild estimates of the number of vertices in each quadrant after CSG operations, which could help rule out some global tree rewrites)

Evaluating all possible pairs would also be a O(n^2) [...] some threshold [...] Its hard to predict the impact (positive or negative) of various optimization attempts without just implementing and measuring them directly.

So... I couldn't resist, I did a crude heuristic that tries to find pair-wise disjoint geometries in the following neigbhours, examining up to K geometries unsuccessfully before calling it a day. But that K allowance is reset if it gathered something in its snowball. Since K is a constant, the search isn't quite quadratic, but I'm too lazy to find out if how much more than linear it actually is.

I ran a few tests and it cuts the N=5 overlap=true example by a sizeable amount (full test results + code): the single loop case is 2x faster than baseline, your nested loops trick is 35% faster than baseline.

The code is a bit crude so I'll work on it it to ensure we fast-union the solids in the right order before I commit anything, but I'm pleased :-D

Btw, added some copyright headers to the new files as per my employer's wishes (I'm doing all this just for fun, though), please let me know if that's an issue.

caching intermediate results to see if it could somehow approximate your manual optimization: it resulted in lots of cache evictions, so not a free pass just yet; left it out of the last batch of commits)

I don't think this type of optimization will be possible, due to the way caching works. Cache keys are basically just the substring / subtree of the CSG Tree which you see from "Design -> Display CSG Tree..."

What I did is to retain a sorted set of the cache keys of the geometries that have been unioned in the LazyGeometry, so that I can create a semantic key ("union(){" + join(componentKeys) + "}) that's order invariant (order of children and/or operations, e.g. same key for union{ a; union { b; c; } } as for union{ a; b; c; } or for union{ c; b; a; }, especially if we flatten the unions with the other PR). An issue though is the priority queue would need to dynamically prioritize already cached subsets of the operands (should be possible?)

I haven't looked enough at the cache size issue but maybe with persistent caching it could become okay to store all that extra data?

keep in mind that CSGTreeNormalizer already basically converts to a sum of products form. You can see CSGTreeNormalizer::match_and_replace

I saw this too early in my discovery of the codebase, concluded it wasn't used in rendering and completely forgot about it! 🤦‍♂️

This pull request introduces 1 alert when merging be8ab70 into 99cae25 - view on LGTM.com

new alerts:

• 1 for Inconsistent definition of copy constructor and assignment ('Rule of Two')

Maybe the LazyGeometry functions can be used to replace this instead? I definitely would prefer more general solutions that can be applied in multiple places over having tons of niche classes for every little thing.

I've seen it, and agree! I've split this logic to a mixed_cache.h file as a starting point

(oh, I think I now see what you mean, did you suggest also moving the on-demand conversion logic to that centralized cache facade? lots more code could get speed ups with that, e.g. intersections / differences)

This pull request introduces 1 alert when merging 7b6e2d2 into 99cae25 - view on LGTM.com

new alerts:

• 1 for Inconsistent definition of copy constructor and assignment ('Rule of Two')

I think could potentially be used beyond the union case
...
(oh, I think I now see what you mean, did you suggest also moving the on-demand conversion logic to that centralized cache facade? lots more code could get speed ups with that, e.g. intersections / differences)

Yeah, I guess I was just meaning that it would be a shame if on demand conversion code ended up not being used beyond the union case, and that it sounds like there's now two places where we have some code choosing between two caches in one way or another. Definitely needs some unified solution, I had at one point looked into making a single GeometryCache which took take a tuple of shared_ptrs for the different types, not sure if that's ideal either, but it was one idea.

Just advocating the DRY principle basically, after all I'm quite proud that my sum contributions to OpenSCAD is in the negative LoC :D

What I did is to retain a sorted set of the cache keys of the geometries that have been unioned in the LazyGeometry, so that I can create a semantic key ("union(){" + join(componentKeys) + "}) that's order invariant (order of children and/or operations, e.g. same key for union{ a; union { b; c; } } as for union{ a; b; c; } or for union{ c; b; a; }, especially if we flatten the unions with the other PR). An issue though is the priority queue would need to dynamically prioritize already cached subsets of the operands (should be possible?)

I see, that's clever and maybe useful in some cases, but I'm not sure. I still don't see it helping with the original script's translate([i,j,0]) ...; because that single loop's translate* never has any repeating values for the combined [i,j,0] ... every single translate is unique. Does that make sense?

* Well its technically a multmatrix node, after being formatted as CSG Tree text.

I haven't looked enough at the cache size issue but maybe with persistent caching it could become okay to store all that extra data?

If there is still some use case for caching intermediate results like that, then maybe instead of appending disjoint geometries to Polysets, we could just useGeometryList (or maybe a new derived DisjointGeometryList?) so that only the shared_ptrs are copied.

As an aside, if you are interested in a more immediate / informal discussion, I'd invite you to join the IRC channel: #openscad on irc.freenode.net
The three of us: t-paul, rcolyer, and I are pretty much on there 24/7 along with some other regular users, and happy to discuss new ideas for improving the project, help with scad scripts, etc.

Ok, so I'm currently reworking the code to be faster, and less complex.

My current approach is to do an analysis phase where I partition union operands into clusters of mutually non-overlapping geometries (no LazyGeometry class anymore). Some of these clusters will be of size 1. Then I reduce each of these non-overlapping clusters by appending in-place to a PolySet (avoids many copies, especially in pathological cases with, say, hundreds of thousands of small geometries). At the end of each cluster reduction I do a single check of integrity and abort that cluster if it's not closed or valid ("explode" it to clusters of size 1 for subsequent normal union), to amortize the costs of integrity checks (which should rarely fail anyway?). The resulting clusters of size 1 are slow-unioned as usual.

I've got a working proof of concept (using CGAL::spatial_sort and CGAL::Union_find) and I'm working on switching to CGAL::box_intersection_d for faster bbox tests.

All this to say, please hold off before reviewing too much of this as I hope the code will become smaller (will also try and merge some of the LazyGeometry methods with existing cache helpers as suggested by @thehans ).

So, that algorithm wasn't such a good idea after all, it did find large patches of fast-unionable polysets and that bit was fast, but if there were any overlaps left it ended up unioning much larger nef polyhedra at the end, making things slower for some large models.

I've now pursued a different route in #3641, which is to use the Polygon Mesh Processing library (and to salvage the fast-union idea from this PR, but that's now almost a detail).

That PR completely obsoletes this one so I'll close this one :-)