from collections import defaultdict
from dataclasses import dataclass, field

import numpy
import networkx as nx

import ifcopenshell
import ifcopenshell.geom
import ifcopenshell.util.element
import ifcopenshell.util.placement
import ifcopenshell.util.unit


cfg = {
    # Mapping of space LongName fragments to category
    "space_mapping": {
        "communication": "technical_space",
        "electrical": "technical_space",
        "equipment": "technical_space",
        "fan": "technical_space",
        "fire": "technical_space",
        "riser": "technical_space",
        "machine room": "technical_space",
        "mechanical": "technical_space",
        "meter room": "technical_space",
        "plant room": "technical_space",
        "services": "technical_space",
        "ventilation": "technical_space",
        "comm.": "technical_space",
        "server": "technical_space",
        "service": "technical_space",
        "shaft": "technical_space",
        "plenum": "technical_space",
        "duct": "technical_space",
        "crac": "technical_space",
        "mep": "technical_space",
        "hvac": "technical_space",
        "room": "habitable_space",
        "interaction": "habitable_space",
        "living": "habitable_space",
        "bed": "habitable_space",
        "kitchen": "habitable_space",
        "office": "habitable_space",
        "reception": "habitable_space",
        "retail": "habitable_space",
        "restaurant": "habitable_space",
        "auditorium": "habitable_space",
        "studio": "habitable_space",
        "external": "outdoor_space",
        "exterior": "outdoor_space",
        "outdoor": "outdoor_space",
        "balcony": "outdoor_space",
        "terrace": "outdoor_space",
        "courtyard": "outdoor_space",
        "patio": "outdoor_space",
        "garden": "outdoor_space",
        "deck": "outdoor_space",
        "roof": "outdoor_space",
        "egress": "circulation_space",
        "escape": "circulation_space",
        "entrance": "circulation_space",
        "corridor": "circulation_space",
        "hallway": "circulation_space",
        "airlock": "circulation_space",
        "foyer": "circulation_space",
        "lobby": "circulation_space",
        "entry": "circulation_space",
        "stair": "circulation_space",
        "lift": "circulation_space",
        "elevator": "circulation_space",
        "circulation": "circulation_space",
        "store": "functional_space",
        "janitor": "functional_space",
        "jan.": "functional_space",
        "cleaner": "functional_space",
        "storage": "functional_space",
        "cupboard": "functional_space",
        "closet": "functional_space",
        "pantry": "functional_space",
        "laundry": "functional_space",
        "utility": "functional_space",
        "bathroom": "bathroom",
        "washroom": "bathroom",
        "restroom": "bathroom",
        "lavatory": "bathroom",
        "shower": "bathroom",
        "toilet": "bathroom",
        "bath": "bathroom",
        "wc": "bathroom",
        "unassigned": "residual_space",
        "unspecified": "residual_space",
        "residual": "residual_space",
        "void": "residual_space"
    },
    # Constituent categories in area types
    "hierarchy": {
        "habitable_area": ["habitable_space"],
        "outdoor_area": ["outdoor_space"],
        "net_floor_area": ["habitable_space", "functional_space", "residual_space", "technical_space", "outdoor_space"],
        "residual_area": ["bathroom", "residual_space", "circulation_space"],
        "usable_floor_area": ["habitable_space", "bathroom", "residual_space", "circulation_space"]
    }
}

LINEAR_TOL = 1.0e-5
WALL_NORMAL_TOL = 1.0e-1
WALL_DEPTH = (1.0e-3, 0.4)
AVG_PLANE_TOL = (1.0e-4, 1.0e-2)

@dataclass(eq=False)
class Face:
    eq: tuple
    exact: object
    pts: numpy.ndarray
    bbox_min: numpy.ndarray
    bbox_max: numpy.ndarray
    shape: int
    part: int
    index: int


@dataclass
class Shape:
    eid: int
    storey: int | None
    parts: list = field(default_factory=list)
    faces: list[Face] = field(default_factory=list)


@dataclass
class Model:
    file: ifcopenshell.file
    shapes: list[Shape] = field(default_factory=list)
    face_count: int = 0
    volume_unit: object = None


@dataclass
class Part:
    shape: object
    faces: list[Face]


def run(src : ifcopenshell.file, category : str) -> ifcopenshell.file:
    model = read_model(src)
    rules = sorted(cfg.get("space_mapping", {}).items(), key=lambda item: len(item[0]), reverse=True)
    wanted_space_kinds = set(cfg["hierarchy"].get(category, ()))
    spaces = set()
    for i, shape in enumerate(model.shapes):
        elem = model.file[shape.eid]
        if not elem.is_a("IfcSpace"):
            continue
        long_name = str(getattr(elem, "LongName", "") or "").casefold()
        kind = None
        for query, mapped_kind in rules:
            if query.casefold() in long_name:
                kind = mapped_kind
                break
        if kind in wanted_space_kinds:
            spaces.add(i)

    walls = selected_walls(model)
    touches = find_touches(model, spaces, walls)
    sides = [side for wall in sorted(walls) for side in wall_sides(model, wall)]

    trim_narrow(model)
    steps = plan_merges(model, spaces, touches, sides)
    graph = nx.Graph()
    graph.add_nodes_from(spaces)
    graph.add_edges_from((a.shape, b.shape) for a, b in steps if a.shape != b.shape)
    components = [sorted(comp) for comp in nx.connected_components(graph)]

    canonicalize_planes(model, [comp for comp in components if len(comp) > 1])

    solids = apply_merges(model, spaces, steps)
    return create_new_ifc(model, components, solids)


def plane_orientation(eq):
    return max(eq, key=abs).to_double() >= 0.


def apply_plane_maps(model, maps_by_shape):
    canonical_by_key = cluster_plane_mappings([m for maps in maps_by_shape.values() for m in maps])
    for shape_i in sorted(maps_by_shape):
        old_by_key = {}
        for old, _new in maps_by_shape[shape_i]:
            old_by_key.setdefault(exact_key(old), old)
        for part in model.shapes[shape_i].parts:
            for plane in (facet.plane_equation() for facet in part.shape.facets()):
                key = exact_key(plane)
                if key in canonical_by_key:
                    old_by_key.setdefault(key, plane)
        resolved = resolved_shape_mappings(old_by_key, canonical_by_key)
        if not resolved:
            continue

        froms, tos = map(list, zip(*resolved))
        for part in model.shapes[shape_i].parts:
            part.shape.map(froms, tos)

    return canonical_by_key


def resolved_shape_mappings(old_by_key, canonical_by_key):
    resolved = {}
    for key, old in sorted(old_by_key.items(), key=lambda item: item[0]):
        canonical = canonical_by_key.get(key)
        if canonical is None:
            continue
        if exact_key(old) != exact_key(canonical):
            resolved.setdefault(key, (old, canonical))
    return list(resolved.values())


def cluster_plane_mappings(mappings):
    parent = {}
    plane_by_key = {}

    def add(eq):
        key = exact_key(eq)
        parent.setdefault(key, key)
        plane_by_key.setdefault(key, eq)
        return key

    def find(key):
        while parent[key] != key:
            parent[key] = parent[parent[key]]
            key = parent[key]
        return key

    def union(a, b):
        ra = find(a)
        rb = find(b)
        if ra == rb:
            return
        parent[max(ra, rb)] = min(ra, rb)

    for old, new in mappings:
        old_key = add(old)
        new_key = add(new)
        if old_key != new_key:
            union(old_key, new_key)

    groups = defaultdict(list)
    for key in sorted(parent):
        groups[find(key)].append(key)

    canonical_by_key = {}
    seen_roots = set()
    for root in sorted(groups, key=lambda r: (len(groups[r]), groups[r][0])):
        if root in seen_roots:
            continue
        roots = {root}
        for key in groups[root]:
            neg_key = exact_key(-plane_by_key[key])
            if neg_key in parent:
                roots.add(find(neg_key))
        seen_roots.update(roots)

        keys = sorted({key for paired_root in roots for key in groups[paired_root]})
        if len(keys) < 2:
            continue
        canonical = canonical_unoriented_cluster_plane([plane_by_key[key] for key in keys])
        for key in keys:
            plane = plane_by_key[key]
            canonical_by_key[key] = canonical if plane_orientation(plane) > 0 else -canonical

    return canonical_by_key


def canonical_unoriented_cluster_plane(planes):
    positives = [plane if plane_orientation(plane) > 0 else -plane for plane in planes]
    units = [unit_plane(plane.to_double()) for plane in positives]
    mean = numpy.average(units, axis=0)

    def distance(item):
        plane, unit = item
        normal_delta = numpy.linalg.norm(unit[:3] - mean[:3]) / AVG_PLANE_TOL[0]
        offset_delta = abs(float(unit[3] - mean[3])) / AVG_PLANE_TOL[1]
        return normal_delta * normal_delta + offset_delta * offset_delta, exact_key(plane)

    return min(zip(positives, units), key=distance)[0]


def exact_key(eq):
    return tuple(v.to_string() for v in eq)


def canonical_unit(eq):
    unit = unit_plane(eq.to_double())
    return tuple(unit if plane_orientation(eq) > 0 else -unit)


def selected_walls(model):
    candidates = []
    for i, shape in enumerate(model.shapes):
        elem = model.file[shape.eid]
        if not elem.is_a("IfcWall"):
            continue

        load_bearing = None
        for pset in ifcopenshell.util.element.get_psets(elem).values():
            if "LoadBearing" in pset:
                load_bearing = pset["LoadBearing"]
                break
        candidates.append((i, load_bearing))

    use_load_property = any(value is not None for _, value in candidates)
    return {i for i, value in candidates if not use_load_property or value is False}


def plan_merges(model, spaces, touches, sides):
    steps = []
    seen = set()

    def add(face_a, face_b):
        key = tuple(sorted(((face_a.shape, face_a.index), (face_b.shape, face_b.index))))
        if key in seen:
            return
        seen.add(key)
        steps.append((face_a, face_b))

    for side_a, side_b in sides:
        left = [face for face in sorted(touches.get(side_a, ()), key=lambda f: f.index) if face.shape in spaces]
        right = [face for face in sorted(touches.get(side_b, ()), key=lambda f: f.index) if face.shape in spaces]
        for face_a in left:
            for face_b in right:
                add(face_a, face_b)

    for a in sorted(spaces):
        for face_a in model.shapes[a].faces:
            for face_b in sorted(touches.get(face_a, ()), key=lambda f: f.index):
                b = face_b.shape
                if b in spaces and a < b:
                    add(face_a, face_b)

    return steps


def canonicalize_planes(model, components):
    maps_by_shape = defaultdict(list)
    face_updates = set()

    for comp in components:
        faces = []
        parent = {}
        units = {}
        buckets = defaultdict(list)

        def find(face_i):
            while parent[face_i] != face_i:
                parent[face_i] = parent[parent[face_i]]
                face_i = parent[face_i]
            return face_i

        def union(a, b):
            ra = find(a)
            rb = find(b)
            if ra != rb:
                parent[ra if ra.index > rb.index else rb] = rb if ra.index > rb.index else ra

        for shape_i in sorted(comp):
            for face in model.shapes[shape_i].faces:
                eq = canonical_unit(face.exact)
                faces.append(face)
                parent[face] = face
                units[face] = eq
                key = plane_bucket(eq, AVG_PLANE_TOL[1], AVG_PLANE_TOL[0])
                for nearby in nearby_buckets(key):
                    for other in buckets.get(nearby, ()):
                        ref = units[other]
                        if (
                            numpy.linalg.norm(numpy.asarray(eq[:3]) - numpy.asarray(ref[:3])) <= AVG_PLANE_TOL[0]
                            and abs(float(eq[3]) - float(ref[3])) <= AVG_PLANE_TOL[1]
                        ):
                            union(face, other)
                buckets[key].append(face)

        clusters = defaultdict(list)
        for face in faces:
            clusters[find(face)].append(face)

        for cluster in clusters.values():
            if len(cluster) < 2:
                continue
            cluster = sorted(cluster, key=lambda f: f.index)

            for face in cluster:
                maps_by_shape[face.shape].append((face.exact, face.exact))
                face_updates.add(face)
            for a, b in zip(cluster, cluster[1:]):
                maps_by_shape[a.shape].append((a.exact, b.exact))

    canonical_by_key = apply_plane_maps(model, maps_by_shape)

    for face in sorted(face_updates, key=lambda f: f.index):
        exact = canonical_by_key.get(exact_key(face.exact))
        if exact is not None:
            face.exact = exact
            face.eq = exact.to_double()

def read_model(file):
    model = Model(file)

    settings = ifcopenshell.geom.settings(
        USE_WORLD_COORDS=True,
        WELD_VERTICES=False,
        DISABLE_OPENING_SUBTRACTIONS=True,
        ITERATOR_OUTPUT=ifcopenshell.ifcopenshell_wrapper.NATIVE
    )
    for item in ifcopenshell.geom.iterate(settings, file, geometry_library="cgal", include=("IfcSpace", "IfcWall")):
        for i in range(item.geometry.size()):
            try:
                parts = list(item.geometry.item(i).convex_decomposition())
            except Exception:
                continue
            if not parts:
                continue

            elem = file[item.id]
            shape = Shape(item.id, shape_storey(elem))
            model.shapes.append(shape)

            shape_i = len(model.shapes) - 1
            for part in parts:
                add_part(model, shape_i, part)
    return model


def add_part(model, shape_i, part):
    shape = model.shapes[shape_i]
    part_i = len(shape.parts)
    faces = part_faces(part, shape_i, part_i)
    for face in faces:
        face.index = model.face_count
        model.face_count += 1
        shape.faces.append(face)
    shape.parts.append(Part(part, faces))


def part_faces(part, shape_i, part_i):
    # Request facets before halfspaces; some wrapper builds expose unstable
    # facet enumeration if the halfspace view is requested first.
    facets = part.facets()
    half = part.halfspaces()
    axes = [f.axis() for f in facets]
    positions = [f.position() for f in facets]
    facet_pts = [tuple(v.position().to_double() for v in f.vertices()) for f in facets]
    facet_keys = [
        tuple(-v for v in axis.to_double()) + (axis.dot(position).to_double(),)
        for axis, position in zip(axes, positions)
    ]

    faces = []
    for exact in (facet.plane_equation() for facet in half.facets()):
        eq = exact.to_double()
        matches = [
            facet_pts[i]
            for i, key in enumerate(facet_keys)
            if numpy.linalg.norm(unit_plane(eq) - unit_plane(key)) < LINEAR_TOL
        ]
        if not matches:
            raise ValueError("Halfspace face without matching polyhedron facet")
        pts = numpy.array(sum(matches, ()))
        faces.append(Face(eq, exact, pts, pts.min(axis=0), pts.max(axis=0), shape_i, part_i, -1))
    return faces


def find_touches(model, spaces, walls):
    touches = defaultdict(set)

    for _storey, wall_faces, space_faces in touch_partitions(model, spaces, walls):
        index = face_plane_index(space_faces)
        for wall_face in wall_faces:
            fa = wall_face
            for space_face in opposite_face_candidates(index, fa.eq):
                fb = space_face
                if not bounds_overlap(fa.bbox_min, fa.bbox_max, fb.bbox_min, fb.bbox_max, LINEAR_TOL):
                    continue
                if not opposite_planes(fa.eq, fb.eq):
                    continue
                if not overlap_2d(fa.pts, fb.pts, fa.eq[:3], fb.eq[:3]):
                    continue
                touches[wall_face].add(space_face)
                touches[space_face].add(wall_face)

    faces_by_storey = defaultdict(list)
    for shape_i in sorted(spaces):
        faces_by_storey[model.shapes[shape_i].storey].extend(model.shapes[shape_i].faces)
    storey_key = lambda item: (item[0] is None, -1 if item[0] is None else item[0])

    for _storey, faces in sorted(faces_by_storey.items(), key=storey_key):
        index = face_plane_index(faces)
        for face_a in faces:
            fa = face_a
            for face_b in opposite_face_candidates(index, fa.eq):
                if face_a.index >= face_b.index:
                    continue
                fb = face_b
                if fa.shape == fb.shape:
                    continue
                if not bounds_overlap(fa.bbox_min, fa.bbox_max, fb.bbox_min, fb.bbox_max, LINEAR_TOL):
                    continue
                if not opposite_planes(fa.eq, fb.eq):
                    continue
                if not overlap_2d(fa.pts, fb.pts, fa.eq[:3], fb.eq[:3]):
                    continue
                touches[face_a].add(face_b)
                touches[face_b].add(face_a)

    return touches


def touch_partitions(model, spaces, walls):
    space_faces_by_storey = defaultdict(list)
    wall_bounds = {}
    for shape_i in sorted(spaces):
        space_faces_by_storey[model.shapes[shape_i].storey].extend(model.shapes[shape_i].faces)
    for shape_i in sorted(walls):
        wall_bounds[shape_i] = faces_bounds(model.shapes[shape_i].faces)

    key = lambda item: (item[0] is None, -1 if item[0] is None else item[0])
    for storey, space_faces in sorted(space_faces_by_storey.items(), key=key):
        storey_bounds = faces_bounds(space_faces)
        wall_faces = []
        for shape_i in sorted(walls):
            shape = model.shapes[shape_i]
            if (
                storey is None
                or shape.storey is None
                or shape.storey == storey
                or bounds_overlap(*wall_bounds[shape_i], *storey_bounds, LINEAR_TOL)
            ):
                wall_faces.extend(shape.faces)
        if wall_faces and space_faces:
            yield storey, wall_faces, space_faces

def face_plane_index(faces):
    index = defaultdict(list)
    for face in faces:
        index[plane_bucket(face.eq, LINEAR_TOL, LINEAR_TOL)].append(face)
    return index


def opposite_face_candidates(index, eq):
    key = plane_bucket(tuple(-v for v in unit_plane(eq)), LINEAR_TOL, LINEAR_TOL)
    seen = set()
    for nearby in nearby_buckets(key):
        for face_i in index.get(nearby, ()):
            if face_i not in seen:
                seen.add(face_i)
                yield face_i


def wall_sides(model, wall):
    faces = model.shapes[wall].faces
    normal_index = defaultdict(list)
    for face in faces:
        normal_index[normal_bucket(face.eq)].append(face)

    candidates = []
    seen = set()
    for a in faces:
        ea = a.eq
        key = normal_bucket(tuple(-v for v in unit_plane(ea)[:3]))
        for nearby in nearby_buckets(key):
            for b in normal_index.get(nearby, ()):
                if a.index >= b.index or (a, b) in seen:
                    continue
                seen.add((a, b))
                eb = b.eq
                if numpy.linalg.norm(numpy.array(ea[:3]) + numpy.array(eb[:3])) >= WALL_NORMAL_TOL:
                    continue
                depth = ea[3] + eb[3]
                if depth > WALL_DEPTH[0]:
                    candidates.append((depth, a, b))

    sides = []
    for depth, a, b in sorted(candidates, key=lambda item: (item[0], item[1].index, item[2].index)):
        if depth > WALL_DEPTH[1]:
            continue
        if overlap_2d(a.pts, b.pts, a.eq[:3], b.eq[:3]):
            sides.append((a, b))
    return sides


def trim_narrow(model):
    maps = []
    drop = set()
    for shape_i, shape in enumerate(model.shapes):
        for part_i, part in enumerate(shape.parts):
            exacts = [face.exact for face in part.faces]
            neg_keys = {}
            for i, eq in enumerate(exacts):
                neg_keys.setdefault(tuple(-v for v in round_eq(eq)), i)
            for exact in exacts:
                idx = neg_keys.get(round_eq(exact))
                if idx is not None:
                    maps.append((exact, -exacts[idx]))
                    maps.append((-exact, exacts[idx]))
                    drop.add((shape_i, part_i))

    canonical_by_key = cluster_plane_mappings(maps)
    old_by_key = {}
    for old, _new in maps:
        old_by_key.setdefault(exact_key(old), old)
    resolved = resolved_shape_mappings(old_by_key, canonical_by_key)
    froms = [old for old, _new in resolved]
    tos = [new for _old, new in resolved]
    for shape_i, shape in enumerate(model.shapes):
        kept = []
        shape.faces = []
        for part_i, part in enumerate(shape.parts):
            if (shape_i, part_i) in drop:
                continue
            half = part.shape.halfspaces()
            if resolved:
                half.map(froms, tos)
            new_part_i = len(kept)
            for face in part.faces:
                face.part = new_part_i
                update_cached_plane(face, canonical_by_key)
                shape.faces.append(face)
            kept.append(Part(half, part.faces))
        shape.parts = kept


def update_cached_plane(face, canonical_by_key):
    exact = canonical_by_key.get(exact_key(face.exact))
    if exact is None:
        return
    face.exact = exact
    face.eq = exact.to_double()


def apply_merges(model, spaces, steps):
    maps_by_shape = defaultdict(list)
    for face_a, face_b in steps:
        eq_a = face_a.exact
        eq_b = face_b.exact
        maps_by_shape[face_a.shape].append((eq_a, -eq_b))
        maps_by_shape[face_b.shape].append((eq_b, -eq_a))

    apply_plane_maps(model, maps_by_shape)
    return {shape_i: shape_solid(model, shape_i) for shape_i in sorted(spaces)}


def shape_solid(model, shape_i):
    parts = model.shapes[shape_i].parts
    if not parts:
        return None
    if len(parts) == 1:
        return parts[0].shape.solid()

    faces_by_part = [part.faces for part in parts]
    graph = shared_part_graph(faces_by_part)
    groups = ordered_part_groups(len(parts), graph)
    solids = [union_solids(parts[part_i].shape.solid() for part_i in group) for group in groups]
    return union_solids(solids)


def shared_part_graph(faces_by_part):
    index = defaultdict(list)
    for part_i, faces in enumerate(faces_by_part):
        for face_i, face in enumerate(faces):
            index[plane_bucket(face.eq, LINEAR_TOL, LINEAR_TOL)].append((part_i, face_i))

    graph = defaultdict(set)
    for part_i, faces in enumerate(faces_by_part):
        for face in faces:
            key = plane_bucket(tuple(-v for v in unit_plane(face.eq)), LINEAR_TOL, LINEAR_TOL)
            for nearby in nearby_buckets(key):
                for other_part_i, other_face_i in index.get(nearby, ()):
                    if part_i >= other_part_i:
                        continue
                    other = faces_by_part[other_part_i][other_face_i]
                    if not bounds_overlap(face.bbox_min, face.bbox_max, other.bbox_min, other.bbox_max, LINEAR_TOL):
                        continue
                    if not opposite_planes(face.eq, other.eq):
                        continue
                    if not overlap_2d(face.pts, other.pts, face.eq[:3], other.eq[:3]):
                        continue
                    graph[part_i].add(other_part_i)
                    graph[other_part_i].add(part_i)
    return graph


def ordered_part_groups(count, graph):
    seen = set()
    groups = []
    for start in range(count):
        if start in seen:
            continue
        seen.add(start)
        order = [start]
        queue = [start]
        for part_i in queue:
            for other in sorted(graph.get(part_i, ())):
                if other in seen:
                    continue
                seen.add(other)
                order.append(other)
                queue.append(other)
        groups.append(order)
    return groups


def create_new_ifc(model, components, solids):
    file = model.file
    made = []
    for comp in components:
        solid, used = component_solid(comp, solids)
        if solid is not None:
            space = make_space(model, used, solid)
            if space is not None:
                made.append(space)

    out = ifcopenshell.file(schema=file.schema)
    out.header.header_data.assign(file.header.header_data)
    for project in file.by_type("IfcProject"):
        out.add(project)
    for rel in file.by_type("IfcRelAggregates"):
        if spatial_node(rel.RelatingObject) and all(spatial_node(obj) for obj in rel.RelatedObjects):
            out.add(rel)

    products = set(made)
    rels = set()
    for product in products:
        out.add(product)
        rels.update(file.get_inverse(product))
    for rel in sorted(rels, key=lambda entity: entity.id()):
        if rel.is_a("IfcRelAggregates"):
            add_filtered_relation(out, rel, "RelatedObjects", products)
        elif rel.is_a("IfcRelContainedInSpatialStructure"):
            add_filtered_relation(out, rel, "RelatedElements", products)
        elif rel.is_a("IfcRelDefinesByProperties"):
            add_filtered_relation(out, rel, "RelatedObjects", products)
    return out


def component_solid(comp, solids):
    used = [shape_i for shape_i in comp if shape_i in solids]
    return union_solids(solids[shape_i] for shape_i in used), used


def union_solids(solids):
    solids = list(solids)
    if not solids:
        return None
    if len(solids) == 1:
        return solids[0]
    return ifcopenshell.ifcopenshell_wrapper.nary_union(solids)


def spatial_node(elem):
    return elem.is_a("IfcProject") or (
        not elem.is_a("IfcSpace")
        and (elem.is_a("IfcSpatialStructureElement") or elem.is_a("IfcSpatialElement"))
    )


def add_filtered_relation(dst, rel, related_attr, products):
    related = list(getattr(rel, related_attr))
    kept = [obj for obj in related if obj in products]
    if not kept:
        return
    if len(kept) == len(related):
        dst.add(rel)
        return

    owner = dst.add(rel.OwnerHistory) if rel.OwnerHistory else None
    if rel.is_a("IfcRelAggregates"):
        dst.createIfcRelAggregates(
            ifcopenshell.guid.new(), owner, rel.Name, rel.Description, dst.add(rel.RelatingObject), [dst.add(o) for o in kept]
        )
    elif rel.is_a("IfcRelContainedInSpatialStructure"):
        dst.createIfcRelContainedInSpatialStructure(
            ifcopenshell.guid.new(), owner, rel.Name, rel.Description, [dst.add(o) for o in kept], dst.add(rel.RelatingStructure)
        )
    elif rel.is_a("IfcRelDefinesByProperties"):
        dst.createIfcRelDefinesByProperties(
            ifcopenshell.guid.new(),
            owner,
            rel.Name,
            rel.Description,
            [dst.add(o) for o in kept],
            dst.add(rel.RelatingPropertyDefinition),
        )


def make_space(model, comp, solid):
    file = model.file
    volume = solid.volume().to_double()
    verts, faces = solid_faces(solid)
    verts = numpy.array(verts).reshape((-1, 3)) / ifcopenshell.util.unit.calculate_unit_scale(file)

    first = next(iter(comp))
    orig = file[model.shapes[first].eid]
    rep = getattr(orig, "Representation", None)
    bodies = [] if rep is None else [r for r in rep.Representations if r.RepresentationIdentifier == "Body"]
    if not bodies:
        return None

    parent = parent_of(orig)
    if parent is not None:
        inv = numpy.linalg.inv(ifcopenshell.util.placement.get_local_placement(parent.ObjectPlacement))
        verts = numpy.concatenate((verts, numpy.ones((len(verts), 1))), axis=1)
        verts = numpy.array([(inv @ p)[:3] for p in verts])

    points = [file.createIfcCartesianPoint(p.tolist()) for p in verts]
    brep_faces = [
        file.createIfcFace([file.createIfcFaceOuterBound(file.createIfcPolyLoop([points[i] for i in face]), True)])
        for face in faces
    ]
    rep = file.createIfcShapeRepresentation(*list(bodies[0])[:2], "Brep", [file.createIfcFacetedBrep(file.createIfcClosedShell(brep_faces))])
    space = file.createIfcSpace(
        ifcopenshell.guid.new(),
        orig.OwnerHistory,
        LongName=" ".join(sorted(filter(None, (getattr(file[model.shapes[i].eid], "LongName", None) for i in comp)))),
        ObjectPlacement=file.createIfcLocalPlacement(
            parent.ObjectPlacement if parent is not None else None,
            file.createIfcAxis2Placement3D(file.createIfcCartesianPoint((0.0, 0.0, 0.0))),
        ),
        Representation=file.createIfcProductDefinitionShape(Representations=[rep]),
    )

    if model.volume_unit is None:
        model.volume_unit = file.createIfcSIUnit(UnitType="VOLUMEUNIT", Name="CUBIC_METRE")
    file.createIfcRelDefinesByProperties(
        ifcopenshell.guid.new(),
        orig.OwnerHistory,
        RelatedObjects=[space],
        RelatingPropertyDefinition=file.createIfcElementQuantity(
            ifcopenshell.guid.new(),
            orig.OwnerHistory,
            Name="Qto_SpaceBaseQuantities",
            Quantities=[file.createIfcQuantityVolume(Name="GrossVolume", Unit=model.volume_unit, VolumeValue=volume)],
        ),
    )
    if parent is not None:
        file.createIfcRelAggregates(ifcopenshell.guid.new(), orig.OwnerHistory, RelatingObject=parent, RelatedObjects=[space])
    return space


def parent_of(elem):
    return elem.Decomposes[0].RelatingObject if getattr(elem, "Decomposes", ()) else None


def shape_storey(elem):
    if elem.is_a("IfcBuildingStorey"):
        return elem.id()
    for candidate in (parent_of(elem), ifcopenshell.util.element.get_container(elem)):
        seen = set()
        while candidate is not None and candidate.id() not in seen:
            seen.add(candidate.id())
            if candidate.is_a("IfcBuildingStorey"):
                return candidate.id()
            candidate = parent_of(candidate)
    return None


def solid_faces(solid):
    verts = []
    faces = []
    for facet in solid.facets():
        points = [vertex.position().to_double() for vertex in facet.vertices()]
        face = []
        for point in points:
            face.append(len(verts))
            verts.append(point)
        if len(face) >= 3:
            faces.append(tuple(face))
    return verts, faces


def faces_bounds(faces):
    return (
        numpy.vstack([face.bbox_min for face in faces]).min(axis=0),
        numpy.vstack([face.bbox_max for face in faces]).max(axis=0),
    )


def bounds_overlap(a_min, a_max, b_min, b_max, tol):
    return bool(numpy.all(a_min <= b_max + tol) and numpy.all(b_min <= a_max + tol))


def plane_bucket(eq, offset_tol, normal_tol):
    eq = unit_plane(eq)
    return tuple(int(round(float(v) / normal_tol)) for v in eq[:3]) + (int(round(float(eq[3]) / offset_tol)),)


def normal_bucket(eq):
    return tuple(int(round(float(v) / WALL_NORMAL_TOL)) for v in unit_plane(eq)[:3])


def nearby_buckets(key):
    for i in range(3 ** len(key)):
        n = i
        step = []
        for _ in key:
            step.append(n % 3 - 1)
            n //= 3
        yield tuple(k + delta for k, delta in zip(key, step))


def opposite_planes(a, b):
    return numpy.linalg.norm(numpy.array(a) + numpy.array(b)) < LINEAR_TOL


def unit_plane(eq):
    eq = numpy.asarray(eq, dtype=float)
    n = numpy.linalg.norm(eq[:3])
    return eq if n == 0.0 else eq / n


def overlap_2d(pa, pb, na, nb):
    for a, b, n in ((pa, pb, na), (pb, pa, nb)):
        a = numpy.asarray(a, dtype=float)
        b = numpy.asarray(b, dtype=float)
        n = unit_plane(tuple(n) + (0.0,))[:3]
        center = numpy.average(a, axis=0)
        local_a = a - center
        local_b = b - center - numpy.outer((b - center) @ n, n)
        axis = numpy.asarray(min(([0, 0, 1], [1, 0, 0]), key=lambda c: abs(numpy.dot(n, c))), dtype=float)
        y = numpy.cross(n, axis)
        y /= numpy.linalg.norm(y)
        x = numpy.cross(y, n)
        a2 = numpy.column_stack((local_a @ x, local_a @ y))
        b2 = numpy.column_stack((local_b @ x, local_b @ y))
        overlap = numpy.minimum(a2.max(axis=0), b2.max(axis=0)) - numpy.maximum(a2.min(axis=0), b2.min(axis=0))
        if not numpy.all(overlap > LINEAR_TOL):
            return False
    return True


def round_eq(eq):
    return tuple(int(round(v / LINEAR_TOL)) for v in eq.to_double())