diff --git a/tools/ci/measure-reference-models.sh b/tools/ci/measure-reference-models.sh index 07caae1a..c51071e0 100755 --- a/tools/ci/measure-reference-models.sh +++ b/tools/ci/measure-reference-models.sh @@ -258,7 +258,7 @@ finished() { # dir name started # longest ones do not start last (the BIM4LCA office's mesh_bodies_vs_ifc took # 477 s of its 1,953 run one after another, run 41589642386). model_tasks() { # name rvt ref - local name="$1" rvt="$2" ref="$3" s + local name="$1" rvt="$2" ref="$3" s meshes= echo "1 base" # ports_vs_ifc reads the findings the schema step's validation wrote, so it # runs after it, in the same task. @@ -270,8 +270,10 @@ model_tasks() { # name rvt ref wanted "$s" || continue [ "$s" = ports_vs_ifc ] && wanted schema && continue case "$s" in - mesh_bodies_vs_ifc) echo "10 ifc:$s" ;; - saved_meshes_vs_ifc) echo "8 ifc:$s" ;; + # These two mesh the same bodies of the same two files, so they run + # one after the other as one task and the second reads what the + # first meshed. + mesh_bodies_vs_ifc | saved_meshes_vs_ifc) meshes=1 ;; layer_sets_vs_ifc | wall_heights_vs_ifc | door_bodies_vs_ifc) echo "2 ifc:$s" ;; *) echo "1 ifc:$s" ;; esac @@ -282,6 +284,7 @@ model_tasks() { # name rvt ref echo "1 ifc:$s" done fi + [ -n "$meshes" ] && echo "18 meshes" [ -f "$MODELS/${rvt%.rvt}.revit-sets.ifc" ] && wanted revit_property_sets_vs_ifc && echo "2 sets" return 0 } @@ -326,6 +329,14 @@ model_task() { # name rvt ref task s=${task#ifc:} score "$dir/$s.txt" "$SCORERS/$s.py" "$dir/model.ifc" "$MODELS/$ref" ;; + meshes) + for s in saved_meshes_vs_ifc mesh_bodies_vs_ifc; do + wanted "$s" || continue + RVT_MESH_CACHE="$dir/mesh-cache" score "$dir/$s.txt" "$SCORERS/$s.py" "$dir/model.ifc" "$MODELS/$ref" + done + echo "$name: mesh cache $(du -sh "$dir/mesh-cache" 2>/dev/null | cut -f1), removed" + rm -rf "$dir/mesh-cache" + ;; sets) # RE-210: where Revit's export with its own property sets is beside the # model (.revit-sets.ifc, the R2025 models), an export with @@ -350,6 +361,7 @@ model_task() { # name rvt ref task model_done() { # name local dir ms dir=$(model_dir "$1") + rm -rf "$dir/mesh-cache" rm -f "$dir/model.glb" "$dir/validation.jsonl" [ -n "${KEEP_IFC:-}" ] || rm -f "$dir/model.ifc" ms=$(awk '{ms += $2} END {print ms + 0}' "$dir/timings.log" 2>/dev/null || echo 0) diff --git a/tools/re/ifc_common.py b/tools/re/ifc_common.py index 62c686e9..06a6aca8 100644 --- a/tools/re/ifc_common.py +++ b/tools/re/ifc_common.py @@ -1,7 +1,15 @@ """Helpers the tools/re scorers share, one copy of each so they cannot drift.""" +import concurrent.futures +import concurrent.futures.process +import hashlib +import multiprocessing +import os import re +import sys +import time +import ifcopenshell.geom import ifcopenshell.util.placement import numpy as np @@ -36,3 +44,134 @@ def report_mesh_error(error): if message not in MESH_ERRORS: MESH_ERRORS.add(message) print(f"IfcOpenShell cannot mesh: {message}") + + +_MESHING = {} +# Instance id -> the first line of IfcOpenShell's reason, for each part +# `mesh_parts` could not mesh (a part kept in the cache keeps its reason too). +# Nothing prints them: a scorer's output does not change when a part fails. +MESH_FAILURES = {} +# Fewer parts than this are meshed in this process: a forked worker parses +# again the parts of the file it meshes, which a handful do not repay. +POOL_FROM = 16 + + +def _in_worker(): + _MESHING["parent"] = os.getppid() + + +def _mesh_one(part_id): + """(vertices, faces) of the part, or the first line of why IfcOpenShell would not mesh it.""" + # An OOM-killed parent leaves its workers with the queued parts: stop at the next one. + if "parent" in _MESHING and os.getppid() != _MESHING["parent"]: + os._exit(1) + try: + shape = ifcopenshell.geom.create_shape(_MESHING["settings"], _MESHING["model"].by_id(part_id)) + except RuntimeError as error: + return (str(error).strip().splitlines() or [""])[0] or type(error).__name__ + return ( + np.array(shape.geometry.verts, dtype=float).reshape(-1, 3), + np.array(shape.geometry.faces, dtype=int).reshape(-1, 3), + ) + + +def _settings_key(settings): + items = [] + for name in sorted(settings.setting_names()): + try: + items.append(f"{name}={settings.get(name)!r}") + except RuntimeError: + items.append(f"{name}=unset") + return ";".join(items) + + +def _mesh_cache(path, settings): + root = os.environ.get("RVT_MESH_CACHE") + if not root or not path: + return None + stat = os.stat(path) + key = f"{os.path.realpath(path)} {stat.st_size} {stat.st_mtime_ns} {_settings_key(settings)}" + return os.path.join(root, hashlib.sha256(key.encode()).hexdigest()[:16]) + + +def _read_meshes(path): + """What `_write_meshes` kept: instance id -> (vertices, faces), or the reason IfcOpenShell gave for a part it could not mesh.""" + out = {} + with np.load(path, allow_pickle=False) as kept: + ids, reasons, nverts, nfaces = kept["ids"], kept["reasons"], kept["nverts"], kept["nfaces"] + verts, faces = kept["verts"], kept["faces"] + vstart = np.concatenate([[0], np.cumsum(nverts)]) + fstart = np.concatenate([[0], np.cumsum(nfaces)]) + for k, part_id in enumerate(ids): + out[int(part_id)] = str(reasons[k]) or (verts[vstart[k] : vstart[k + 1]], faces[fstart[k] : fstart[k + 1]]) + return out + + +def _write_meshes(path, meshes): + ids = list(meshes) + parts = [(np.empty((0, 3)), np.empty((0, 3), dtype=int)) if isinstance(meshes[i], str) else meshes[i] for i in ids] + scratch = f"{path}.tmp" + with open(scratch, "wb") as f: + np.savez( + f, + ids=np.array(ids, dtype=np.int64), + reasons=np.array([meshes[i] if isinstance(meshes[i], str) else "" for i in ids]), + nverts=np.array([len(v) for v, _ in parts], dtype=np.int64), + nfaces=np.array([len(t) for _, t in parts], dtype=np.int64), + verts=np.concatenate([v for v, _ in parts]), + faces=np.concatenate([t for _, t in parts]), + ) + os.replace(scratch, path) + + +def mesh_parts(model, settings, parts, path=None): + """Instance id -> (vertices (N, 3), faces (M, 3)) for each of `parts` that + IfcOpenShell meshes, one `create_shape` per part spread over the CPUs. The + triangles are those `create_shape` gives, in the same order (the + iterator's are the same, but it took twice the CPU time per element in + a measurement on the BIM4LCA office, and no less wall time with 12 threads + than with 4); a part it cannot mesh is left out, with its reason in + `MESH_FAILURES`. + + With `RVT_MESH_CACHE` set to a directory and the model's `path`, what is + meshed is kept there for the next scorer to read, which Measure sets for + the scorers it runs one after another on one pair of files and deletes + afterwards. Each call keeps its own shard, named by process and time, and + the cache is keyed by the file's path, size and modification time and by + every setting.""" + ids = list(dict.fromkeys(part.id() for part in parts)) + cache = _mesh_cache(path, settings) + known = {} + if cache and os.path.isdir(cache): + for name in sorted(os.listdir(cache)): + if name.endswith(".npz"): + known.update(_read_meshes(os.path.join(cache, name))) + todo = [i for i in ids if i not in known] + cpus = len(os.sched_getaffinity(0)) if hasattr(os, "sched_getaffinity") else os.cpu_count() or 1 + workers = min(int(os.environ.get("RVT_MESH_WORKERS") or cpus), len(todo)) + if len(todo) < POOL_FROM: + workers = 1 + _MESHING.update(model=model, settings=settings) + meshes = None + if workers > 1: + # A forked worker flushes the buffered output it inherits when it ends. + sys.stdout.flush() + sys.stderr.flush() + try: + with concurrent.futures.ProcessPoolExecutor( + workers, mp_context=multiprocessing.get_context("fork"), initializer=_in_worker + ) as pool: + meshes = list(pool.map(_mesh_one, todo, chunksize=1)) + except concurrent.futures.process.BrokenProcessPool: + print("mesh_parts: a worker died, so the parts are meshed one after another", file=sys.stderr) + meshes = None + if meshes is None: + meshes = [_mesh_one(i) for i in todo] + _MESHING.clear() + fresh = dict(zip(todo, meshes)) + if cache and fresh: + os.makedirs(cache, exist_ok=True) + _write_meshes(os.path.join(cache, f"{os.getpid()}-{time.time_ns()}.npz"), fresh) + known.update(fresh) + MESH_FAILURES.update({i: known[i] for i in ids if isinstance(known[i], str)}) + return {i: known[i] for i in ids if not isinstance(known[i], str)} diff --git a/tools/re/mesh_bodies_vs_ifc.py b/tools/re/mesh_bodies_vs_ifc.py index eeaffca4..5fcece3d 100644 --- a/tools/re/mesh_bodies_vs_ifc.py +++ b/tools/re/mesh_bodies_vs_ifc.py @@ -36,12 +36,17 @@ import numpy as np sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) -from opening_boxes_vs_ifc import site_frame # noqa: E402 +from ifc_common import mesh_parts, site_frame # noqa: E402 from roof_bodies_vs_ifc import shapes_of # noqa: E402 from saved_meshes_vs_ifc import body_sources # noqa: E402 SOURCE = "saved_graphics_mesh" PAIR_BUDGET = 1_000_000 +# surface_distance measures every pair up to this many, and beyond it prunes. +EXHAUSTIVE = 200_000 +BLOCK = 8 +SEEDS = 4 +SLIVER = 1e-2 # Most elements scored per class (an even stride through them by Tag), and # most vertices of a body tested against the other's surface, so a file # with thousands of mullions is scored in the time Measure allows. @@ -49,13 +54,14 @@ VERTEX_CAP = 400 -def meshes(path, tags=None, source=None): +def meshes(path, tags=None, source=None, model=None): """Tag -> (class, ObjectType, triangles (T, 3, 3) in internal feet).""" - model = ifcopenshell.open(path) + if model is None: + model = ifcopenshell.open(path) frame = site_frame(model) settings = ifcopenshell.geom.settings() settings.set("use-world-coords", True) - out = {} + chosen = [] for element in model.by_type("IfcElement"): if element.is_a("IfcOpeningElement") or not element.Tag: continue @@ -63,19 +69,20 @@ def meshes(path, tags=None, source=None): continue if source is not None and body_sources(element) != source: continue - parts = [] - for part in shapes_of(element): - try: - shape = ifcopenshell.geom.create_shape(settings, part) - except RuntimeError: + chosen.append((element, shapes_of(element))) + shapes = mesh_parts(model, settings, [part for _, parts in chosen for part in parts], path) + out = {} + for element, parts in chosen: + triangles = [] + for part in parts: + if part.id() not in shapes: continue - verts = np.array(shape.geometry.verts, dtype=float).reshape(-1, 3) - faces = np.array(shape.geometry.faces, dtype=int).reshape(-1, 3) + verts, faces = shapes[part.id()] if len(verts) and len(faces): - parts.append(frame(verts)[faces]) - if not parts: + triangles.append(frame(verts)[faces]) + if not triangles: continue - triangles = np.concatenate(parts) + triangles = np.concatenate(triangles) if element.Tag in out: triangles = np.concatenate([out[element.Tag][2], triangles]) out[element.Tag] = (element.is_a(), element.ObjectType or "", triangles) @@ -119,42 +126,148 @@ def find(k): return sorted(shells.values(), key=lambda v: -abs(v)) +def dot(x, y): + return np.einsum("...k,...k->...", x, y) + + def segment_distances(p, a, b): ab = b - a - t = np.einsum("ntk,tk->nt", p - a, ab) / np.maximum(np.einsum("tk,tk->t", ab, ab), 1e-24) + t = dot(p - a, ab) / np.maximum(dot(ab, ab), 1e-24) q = a + np.clip(t, 0.0, 1.0)[..., None] * ab return np.linalg.norm(p - q, axis=-1) -def surface_distance(points, triangles): - """The largest distance from any of `points` to the nearest triangle.""" +def triangle_terms(triangles): a, b, c = triangles[:, 0], triangles[:, 1], triangles[:, 2] ab, ac = b - a, c - a normal = np.cross(ab, ac) length = np.linalg.norm(normal, axis=1) unit = normal / np.maximum(length, 1e-24)[:, None] - d00 = np.einsum("tk,tk->t", ab, ab) - d01 = np.einsum("tk,tk->t", ab, ac) - d11 = np.einsum("tk,tk->t", ac, ac) + d00, d01, d11 = dot(ab, ab), dot(ab, ac), dot(ac, ac) denom = np.where(length > 1e-12, d00 * d11 - d01 * d01, 1.0) - worst = 0.0 - step = max(1, PAIR_BUDGET // len(triangles)) - for start in range(0, len(points), step): - p = points[start : start + step][:, None, :] - ap = p - a - height = np.einsum("ntk,tk->nt", ap, unit) - q = ap - height[..., None] * unit - d20 = np.einsum("ntk,tk->nt", q, ab) - d21 = np.einsum("ntk,tk->nt", q, ac) - v = (d11 * d20 - d01 * d21) / denom - w = (d00 * d21 - d01 * d20) / denom - inside = (v >= -1e-9) & (w >= -1e-9) & (v + w <= 1 + 1e-9) & (length > 1e-12) - plane = np.where(inside, np.abs(height), np.inf) - nearest = np.minimum.reduce( - [plane, segment_distances(p, a, b), segment_distances(p, b, c), segment_distances(p, c, a)] - ).min(axis=1) - worst = max(worst, float(nearest.max())) - return worst + return a, b, c, unit, length, d00, d01, d11, denom + + +def distances(p, terms): + """Distance from points `p` (..., 1, 3) to triangles whose terms broadcast with them.""" + a, b, c, unit, length, d00, d01, d11, denom = terms + ab, ac = b - a, c - a + ap = p - a + height = dot(ap, unit) + q = ap - height[..., None] * unit + d20 = dot(q, ab) + d21 = dot(q, ac) + v = (d11 * d20 - d01 * d21) / denom + w = (d00 * d21 - d01 * d20) / denom + inside = (v >= -1e-9) & (w >= -1e-9) & (v + w <= 1 + 1e-9) & (length > 1e-12) + plane = np.where(inside, np.abs(height), np.inf) + return np.minimum.reduce([plane, segment_distances(p, a, b), segment_distances(p, b, c), segment_distances(p, c, a)]) + + +def morton(centres): + """The order of points along a Z-order curve, so triangles near each other are near in the order.""" + lo, hi = centres.min(axis=0), centres.max(axis=0) + cells = ((centres - lo) / np.maximum(hi - lo, 1e-12) * 1023).astype(np.uint64) + code = np.zeros(len(centres), dtype=np.uint64) + for bit in range(10): + for axis in range(3): + code |= ((cells[:, axis] >> np.uint64(bit)) & np.uint64(1)) << np.uint64(3 * bit + axis) + return np.argsort(code, kind="stable") + + +def surface_distance(points, triangles): + """The largest distance from any of `points` to the nearest triangle. + + Exact, not approximate: for each point the nearest distance is the minimum of + `distances` over a set of triangles that holds every one that could be nearest, + and each of them is measured by the very arithmetic that measures all of them in + a small job (`EXHAUSTIVE` pairs or fewer, and any input that is not finite), so + the result is the same float. + + A triangle is left out for a point only when a lower bound on its distance + exceeds the nearest distance already found for that point plus a slack. The + slack is 1e-6 ft and 4e-9 of the extent, the largest of any coordinate and any + triangle edge. It has to cover two things. `inside` takes a point up to 1e-9 of + two edges outside the triangle in its plane (or 1e-9 past the third) as inside + and measures it by its plane alone, so the distance computed for it falls short + of the distance to the triangle, which the bounds are on, by up to 1e-9 times + the sum of two edges: under 2e-9 of the extent. And the rounding of these + distances, which is many orders of magnitude smaller. + + The bounds hold for every triangle but the slivers (the sine of the angle between + two edges under SLIVER, whose computed plane is not reliable) and the zero-area + ones (which have no plane, so only their box bounds them): the distance to a + block of eight neighbouring triangles' box, then to a triangle's box, and for a + triangle with a plane the distance to its plane combined with that to the side + planes of its edges. Slivers are measured against every point. + """ + terms = triangle_terms(triangles) + n = len(points) + if n * len(triangles) <= EXHAUSTIVE or not (np.isfinite(points).all() and np.isfinite(triangles).all()): + worst = 0.0 + step = max(1, PAIR_BUDGET // len(triangles)) + for start in range(0, n, step): + nearest = distances(points[start : start + step][:, None, :], terms).min(axis=1) + worst = max(worst, float(nearest.max())) + return worst + a, b, c, unit, length, d00, d01, d11, denom = terms + sliver = (length > 1e-12) & (length < SLIVER * np.sqrt(d00 * d11)) + flat = (length > 1e-12) & ~sliver + lo = np.minimum(np.minimum(a, b), c) + hi = np.maximum(np.maximum(a, b), c) + extent = max(float(np.abs(triangles).max()), float(np.abs(points).max()), float(np.sqrt(max(d00.max(), d11.max())))) + slack = 1e-6 + 4e-9 * extent + sides = [] + for start, end in ((a, b), (b, c), (c, a)): + out = np.cross(end - start, unit) + sides.append((start, out / np.maximum(np.linalg.norm(out, axis=1), 1e-24)[:, None])) + order = morton((lo + hi) / 2) + order = np.concatenate([order, np.full(-len(order) % BLOCK, order[-1])]) + blocks = order.reshape(-1, BLOCK) + block_lo, block_hi = lo[blocks].min(axis=1), hi[blocks].max(axis=1) + seeds = min(SEEDS, len(blocks)) + nearest = np.empty(n) + step = max(1, PAIR_BUDGET // len(blocks)) + for start in range(0, n, step): + p = points[start : start + step] + gap = np.maximum(np.maximum(block_lo[None] - p[:, None], p[:, None] - block_hi[None]), 0.0) + bound = np.sqrt(dot(gap, gap)) + rows = np.arange(len(p)) + seed = np.argpartition(bound, seeds - 1, axis=1)[:, :seeds] if seeds < len(blocks) else np.tile(np.arange(seeds), (len(p), 1)) + best = np.full(len(p), np.inf) + for column in range(seeds): + best = np.minimum(best, distances(p[:, None, :], tuple(x[blocks[seed[:, column]]] for x in terms)).min(axis=1)) + bound[rows, seed[:, column]] = np.inf + pair_rows, pair_blocks = np.nonzero(bound <= (best + slack)[:, None]) + done = 0 + while done < len(pair_rows): + stop = int(np.searchsorted(pair_rows, pair_rows[min(len(pair_rows) - 1, done + PAIR_BUDGET // BLOCK)], side="right")) + row = np.repeat(pair_rows[done:stop], BLOCK) + tri = blocks[pair_blocks[done:stop]].reshape(-1) + done = stop + q = p[row] + gap = np.maximum(np.maximum(lo[tri] - q, q - hi[tri]), 0.0) + floor = np.sqrt(dot(gap, gap)) + height = dot(q - a[tri], unit[tri]) + outside = np.zeros(len(tri)) + for corner, side in sides: + outside = np.maximum(outside, dot(q - corner[tri], side[tri])) + floor = np.where(flat[tri], np.maximum(floor, np.sqrt(height * height + outside * outside)), floor) + first = np.flatnonzero(np.r_[True, row[1:] != row[:-1]]) + lowest = np.minimum.reduceat(floor, first) + group = np.repeat(np.arange(len(first)), np.diff(np.r_[first, len(row)])) + pick = np.flatnonzero(floor == lowest[group]) + pick = pick[np.unique(group[pick], return_index=True)[1]] + best[row[pick]] = np.minimum(best[row[pick]], distances(q[pick][:, None, :], tuple(x[tri[pick]][:, None] for x in terms))[:, 0]) + keep = floor <= best[row] + slack + row, tri = row[keep], tri[keep] + np.minimum.at(best, row, distances(p[row][:, None, :], tuple(x[tri][:, None] for x in terms))[:, 0]) + nearest[start : start + step] = best + slivers = np.flatnonzero(sliver) + for start in range(0, len(slivers), BLOCK): + part = slivers[start : start + BLOCK] + nearest = np.minimum(nearest, distances(points[:, None, :], tuple(x[part] for x in terms)).min(axis=1)) + return float(nearest.max()) def sample(triangles): @@ -179,8 +292,8 @@ def main(): chosen.update(tags[::step]) if step > 1: print(f"{cls}: {len(tags)} bodies, every {step}th scored") + ours = meshes(sys.argv[1], tags=chosen, source=SOURCE, model=model) del model - ours = meshes(sys.argv[1], tags=chosen, source=SOURCE) ref = meshes(sys.argv[2], tags=set(ours)) print(f"{SOURCE} bodies scored: {len(ours)}; in Revit's export {len(ref)}") rows = collections.defaultdict(lambda: collections.defaultdict(collections.Counter)) diff --git a/tools/re/saved_meshes_vs_ifc.py b/tools/re/saved_meshes_vs_ifc.py index b20193b1..f5ba01d4 100644 --- a/tools/re/saved_meshes_vs_ifc.py +++ b/tools/re/saved_meshes_vs_ifc.py @@ -35,14 +35,16 @@ import argparse import collections import math -import multiprocessing +import os import statistics +import sys import ifcopenshell import ifcopenshell.geom import numpy as np -from ifc_common import site_frame +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +from ifc_common import mesh_parts, site_frame # noqa: E402 FT = 0.3048 @@ -80,12 +82,14 @@ def world_boxes(path, tags=None, transform=None, wanted=None, frame=True): if not elements: return out to_internal = site_frame(model) if frame else None - it = ifcopenshell.geom.iterator(settings, model, multiprocessing.cpu_count(), include=elements) - if not it.initialize(): - return out - while True: - shape = it.get() - verts = np.array(shape.geometry.verts, dtype=float).reshape(-1, 3) + # The order the geometry iterator returned them in, one thread, which the + # tables' ties (equal counts) follow. + elements.sort(key=lambda e: e.id()) + shapes = mesh_parts(model, settings, elements, path) + for element in elements: + if element.id() not in shapes: + continue + verts = shapes[element.id()][0] if len(verts): if to_internal is not None: points = to_internal(verts) @@ -94,10 +98,7 @@ def world_boxes(path, tags=None, transform=None, wanted=None, frame=True): if transform: points = np.array([transform(p) for p in points]) box = [float(v) for v in points.min(axis=0)] + [float(v) for v in points.max(axis=0)] - element = model.by_id(shape.id) out[element.Tag] = (element.is_a(), box, sources.get(element.Tag)) - if not it.next(): - break return out