formulas.py

1import numpy as np2from fractions import Fraction3from math import comb, ceil, floor, log4from typing import Callable, Dict, Iterable, Tuple5from .errors import MrlyError6from . import binary78# GRID910def grid_lines(number: int, level: int) -> int:11    return number**level1213def grid_squares(number: int, level: int) -> int:14    return (number**2)**level1516def grid_cubes(number: int, level: int) -> int:17    return (number**3)**level1819# CARPET2021def carpet_fill_squares(number: int, level: int) -> int:22    return (number**2-floor(number/2)**2)**level2324def carpet_void_squares(number: int, level: int) -> int:25    return grid_squares(number, level) - carpet_fill_squares(number, level)2627def carpet_fill_cubes(number: int, level: int) -> int:28    E = ceil(number/2)29    O = floor(number/2)30    return (E**3 + 3*O*E**2)**level3132def carpet_void_cubes(number: int, level: int) -> int:33    return grid_cubes(number, level) - carpet_fill_cubes(number, level)3435# NET3637def net_fill_squares(number: int, level: int) -> int:38    return (number**2-floor((number+1)/2)**2)**level3940def net_void_squares(number: int, level: int) -> int:41    return grid_squares(number, level) - net_fill_squares(number, level)4243def net_fill_cubes(number: int, level: int) -> int:44    E = ceil(number/2)45    O = floor(number/2)46    return (O**3 + 3*E*O**2)**level4748def net_void_cubes(number: int, level: int) -> int:49    return grid_cubes(number, level) - net_fill_cubes(number, level)5051# TREE5253def tree_fill_squares(number: int, level: int) -> int:54    return (number*floor((number+1)/2))**level5556def tree_void_squares(number: int, level: int) -> int:57    return grid_squares(number, level) - tree_fill_squares(number, level)5859def tree_fill_cubes(number: int, level: int) -> int:60    return (number * ceil(number/2)**2)**level6162def tree_void_cubes(number: int, level: int) -> int:63    return grid_cubes(number, level) - tree_fill_cubes(number, level)6465# VOID6667def void_fill_squares(number: int, level: int) -> int:68    return (ceil(number**2/2))**level6970def void_void_squares(number: int, level: int) -> int:71    return grid_squares(number, level) - void_fill_squares(number, level)7273def void_fill_cubes(number: int, level: int) -> int:74    E = ceil(number/2)75    O = floor(number/2)76    return (E**3 + O**3)**level7778def void_void_cubes(number: int, level: int) -> int:79    return grid_cubes(number, level) - void_fill_cubes(number, level)8081# LEVEL-SET (THE SYMMETRIC GENUS)8283def _level_set_base(number: int, dimension: int, levels: Iterable[int]) -> int:84    even = ceil(number / 2)85    odd = floor(number / 2)86    chosen = set(int(j) for j in levels)87    return sum(comb(dimension, j) * (odd ** j) * (even ** (dimension - j))88               for j in chosen if 0 <= j <= dimension)8990def level_set_fill(number: int, level: int, dimension: int, levels: Iterable[int]) -> int:91    return _level_set_base(number, dimension, levels) ** level9293def level_set_fill_squares(number: int, level: int, levels: Iterable[int]) -> int:94    return level_set_fill(number, level, 2, levels)9596def level_set_fill_cubes(number: int, level: int, levels: Iterable[int]) -> int:97    return level_set_fill(number, level, 3, levels)9899def level_set_void_squares(number: int, level: int, levels: Iterable[int]) -> int:100    return grid_squares(number, level) - level_set_fill_squares(number, level, levels)101102def level_set_void_cubes(number: int, level: int, levels: Iterable[int]) -> int:103    return grid_cubes(number, level) - level_set_fill_cubes(number, level, levels)104105# RATIOS106107def carpet_ratio_2d(number: int, level: int) -> float:108    return carpet_fill_squares(number, level) / grid_squares(number, level)109110def carpet_ratio_3d(number: int, level: int) -> float:111    return carpet_fill_cubes(number, level) / grid_cubes(number, level)112113def net_ratio_2d(number: int, level: int) -> float:114    return net_fill_squares(number, level) / grid_squares(number, level)115116def net_ratio_3d(number: int, level: int) -> float:117    return net_fill_cubes(number, level) / grid_cubes(number, level)118119def tree_ratio_2d(number: int, level: int) -> float:120    return tree_fill_squares(number, level) / grid_squares(number, level)121122def tree_ratio_3d(number: int, level: int) -> float:123    return tree_fill_cubes(number, level) / grid_cubes(number, level)124125def void_ratio_2d(number: int, level: int) -> float:126    return void_fill_squares(number, level) / grid_squares(number, level)127128def void_ratio_3d(number: int, level: int) -> float:129    return void_fill_cubes(number, level) / grid_cubes(number, level)130131def carpet_ratio_6d(number: int, level: int) -> float:132    return carpet_fill_triangles(number, level) / grid_triangles(number, level)133134def net_ratio_6d(number: int, level: int) -> float:135    return net_fill_triangles(number, level) / grid_triangles(number, level)136137def tree_ratio_6d(number: int, level: int) -> float:138    return tree_fill_triangles(number, level) / grid_triangles(number, level)139140def void_ratio_6d(number: int, level: int) -> float:141    return void_fill_triangles(number, level) / grid_triangles(number, level)142143# DIMENSIONS144145def calculate_dimension(design: Callable, number: int, dimension: int) -> float:146    if number == 1:147        return float(dimension)148    grid = design(number)149    fill = np.sum(grid)150    if fill <= 0:151        return 0.0152    return log(fill) / log(number)153154def carpet_2d_dimension(number: int) -> float:155    return calculate_dimension(binary.carpet_2d, number, 2)156157def carpet_3d_dimension(number: int) -> float:158    return calculate_dimension(binary.carpet_3d, number, 3)159160def net_2d_dimension(number: int) -> float:161    return calculate_dimension(binary.net_2d, number, 2)162163def net_3d_dimension(number: int) -> float:164    return calculate_dimension(binary.net_3d, number, 3)165166def tree_2d_dimension(number: int) -> float:167    return calculate_dimension(binary.tree_2d, number, 2)168169def tree_3d_dimension(number: int) -> float:170    return calculate_dimension(binary.tree_3d, number, 3)171172def void_2d_dimension(number: int) -> float:173    return calculate_dimension(binary.void_2d, number, 2)174175def void_3d_dimension(number: int) -> float:176    return calculate_dimension(binary.void_3d, number, 3)177178# TRIANGLES179180_fill_diag_cache: Dict[str, np.ndarray] = {}181_total_diag_cache: Dict[int, np.ndarray] = {}182183def _fill_by_diag(grid: np.ndarray, n: int) -> np.ndarray:184    max_d = 3 * (n - 1)185    f = np.zeros(max_d + 1, dtype=np.int64)186    for cx in range(n):187        for cy in range(n):188            for cz in range(n):189                if grid[cx, cy, cz]:190                    f[cx + cy + cz] += 1191    return f192193def _total_by_diag(n: int) -> np.ndarray:194    if n in _total_diag_cache:195        return _total_diag_cache[n]196    max_d = 3 * (n - 1)197    t = np.zeros(max_d + 1, dtype=np.int64)198    for cx in range(n):199        for cy in range(n):200            for cz in range(n):201                t[cx + cy + cz] += 1202    _total_diag_cache[n] = t203    return t204205def _diag_convolve(f_prev: np.ndarray, f_base: np.ndarray, n: int) -> np.ndarray:206    len_prev = len(f_prev)207    len_base = len(f_base)208    max_s = n * (len_prev - 1) + (len_base - 1)209    f_new = np.zeros(max_s + 1, dtype=np.int64)210    for a in range(len_prev):211        if f_prev[a] == 0:212            continue213        for b in range(len_base):214            if f_base[b] == 0:215                continue216            f_new[n * a + b] += f_prev[a] * f_base[b]217    return f_new218219def _visited_diags(m: int) -> Tuple[list, list]:220    k = (3 * (4 * m - 1)) // 2221    d = k // 4222    if m % 2 == 0:223        return [d - 1, d], [4, 4]224    return [d - 2, d - 1, d], [1, 6, 1]225226def _cut_counts(design_fn: Callable, number: int, level: int) -> Tuple[int, int]:227    key = f"{design_fn.__name__}_{number}"228    if key not in _fill_diag_cache:229        _fill_diag_cache[key] = _fill_by_diag(design_fn(number), number)230    f_base = _fill_diag_cache[key]231    t_base = _total_by_diag(number)232    f_cur = f_base.copy()233    t_cur = t_base.copy()234    for _ in range(1, level):235        f_cur = _diag_convolve(f_cur, f_base, number)236        t_cur = _diag_convolve(t_cur, t_base, number)237    m = number ** level238    diags, weights = _visited_diags(m)239    fill = 0240    total = 0241    for d, w in zip(diags, weights):242        if d < len(f_cur):243            fill += w * int(f_cur[d])244        if d < len(t_cur):245            total += w * int(t_cur[d])246    return int(fill), int(total - fill)247248def grid_triangles(number: int, level: int) -> int:249    return 6 * number ** (2 * level)250251def carpet_fill_triangles(number: int, level: int) -> int:252    fill, _ = _cut_counts(binary.carpet_3d, number, level)253    return fill254255def carpet_void_triangles(number: int, level: int) -> int:256    _, void = _cut_counts(binary.carpet_3d, number, level)257    return void258259def net_fill_triangles(number: int, level: int) -> int:260    fill, _ = _cut_counts(binary.net_3d, number, level)261    return fill262263def net_void_triangles(number: int, level: int) -> int:264    _, void = _cut_counts(binary.net_3d, number, level)265    return void266267def tree_fill_triangles(number: int, level: int) -> int:268    fill, _ = _cut_counts(binary.tree_3d, number, level)269    return fill270271def tree_void_triangles(number: int, level: int) -> int:272    _, void = _cut_counts(binary.tree_3d, number, level)273    return void274275def void_fill_triangles(number: int, level: int) -> int:276    fill, _ = _cut_counts(binary.void_3d, number, level)277    return fill278279def void_void_triangles(number: int, level: int) -> int:280    _, void = _cut_counts(binary.void_3d, number, level)281    return void282283# SURFACE AREA284285def _visible_faces(occ: np.ndarray) -> int:286    total = 0287    for axis in range(3):288        lo = np.take(np.pad(occ, [(1, 0) if i == axis else (0, 0) for i in range(3)]),289                     range(occ.shape[axis]), axis=axis)290        hi = np.take(np.pad(occ, [(0, 1) if i == axis else (0, 0) for i in range(3)]),291                     range(1, occ.shape[axis] + 1), axis=axis)292        total += int(np.clip(occ - lo, 0, None).sum())293        total += int(np.clip(occ - hi, 0, None).sum())294    return total295296def _vh_state(grid: np.ndarray) -> Tuple[int, int]:297    occ = (grid != 0).astype(np.int64)298    vf = _visible_faces(occ)299    return vf, 6 * int(occ.sum()) - vf300301def _solve_2x2(a, b, c, d, p, q):302    det = a * d - b * c303    x = Fraction(p * d - b * q, det)304    y = Fraction(a * q - p * c, det)305    return x, y306307def _surface_transfer(design_fn: Callable, number: int):308    base = design_fn(number)309    g1 = base310    g2 = np.kron(g1, base).astype(np.uint8)311    g3 = np.kron(g2, base).astype(np.uint8)312    (v0, h0), (v1, h1), (v2, h2) = _vh_state(g1), _vh_state(g2), _vh_state(g3)313    m00, m01 = _solve_2x2(v0, h0, v1, h1, v1, v2)314    m10, m11 = _solve_2x2(v0, h0, v1, h1, h1, h2)315    return ((m00, m01), (m10, m11)), (v0, h0)316317def _surface_general(design_fn: Callable, number: int, level: int) -> int:318    grid = design_fn(number)319    s1, h1 = _vh_state(grid)320    if level == 1:321        return s1322    if h1 == 0:323        return s1 * int((grid != 0).sum()) ** (level - 1)324    M, v1 = _surface_transfer(design_fn, number)325    R = [[Fraction(1), Fraction(0)], [Fraction(0), Fraction(1)]]326    for _ in range(level - 1):327        R = [[R[0][0] * M[0][0] + R[0][1] * M[1][0], R[0][0] * M[0][1] + R[0][1] * M[1][1]],328             [R[1][0] * M[0][0] + R[1][1] * M[1][0], R[1][0] * M[0][1] + R[1][1] * M[1][1]]]329    return int(R[0][0] * v1[0] + R[0][1] * v1[1])330331def carpet_surface(number: int, level: int) -> int:332    return _surface_general(binary.carpet_3d, number, level)333334def net_surface(number: int, level: int) -> int:335    return _surface_general(binary.net_3d, number, level)336337def tree_surface(number: int, level: int) -> int:338    return _surface_general(binary.tree_3d, number, level)339340def void_surface(number: int, level: int) -> int:341    return _surface_general(binary.void_3d, number, level)342343def surface_visible_faces(design_fn: Callable, number: int, level: int) -> int:344    grid = design_fn(number)345    if level > 1:346        base = grid347        for _ in range(1, level):348            grid = np.kron(grid, base).astype(np.uint8)349    occ = (grid != 0).astype(np.int64)350    total = 0351    for axis in range(3):352        lo = np.take(np.pad(occ, [(1, 0) if i == axis else (0, 0) for i in range(3)]),353                     range(occ.shape[axis]), axis=axis)354        hi = np.take(np.pad(occ, [(0, 1) if i == axis else (0, 0) for i in range(3)]),355                     range(1, occ.shape[axis] + 1), axis=axis)356        total += int(np.clip(occ - lo, 0, None).sum())357        total += int(np.clip(occ - hi, 0, None).sum())358    return total359360_GRID_FN = {2: grid_squares, 3: grid_cubes}361362_FILL_FN = {363    (binary.carpet_2d.__name__, 2): carpet_fill_squares,364    (binary.net_2d.__name__, 2): net_fill_squares,365    (binary.tree_2d.__name__, 2): tree_fill_squares,366    (binary.htree_2d.__name__, 2): tree_fill_squares,367    (binary.vtree_2d.__name__, 2): tree_fill_squares,368    (binary.void_2d.__name__, 2): void_fill_squares,369    (binary.carpet_3d.__name__, 3): carpet_fill_cubes,370    (binary.net_3d.__name__, 3): net_fill_cubes,371    (binary.tree_3d.__name__, 3): tree_fill_cubes,372    (binary.xtree_3d.__name__, 3): tree_fill_cubes,373    (binary.ytree_3d.__name__, 3): tree_fill_cubes,374    (binary.ztree_3d.__name__, 3): tree_fill_cubes,375    (binary.void_3d.__name__, 3): void_fill_cubes,376}377378def _shift_occ(occ: np.ndarray, axis: int) -> np.ndarray:379    s = np.zeros_like(occ)380    src = [slice(None)] * occ.ndim381    dst = [slice(None)] * occ.ndim382    src[axis] = slice(1, None)383    dst[axis] = slice(0, -1)384    s[tuple(dst)] = occ[tuple(src)]385    return s386387def _core_edges_grid(grid: np.ndarray) -> int:388    occ = (grid != 0)389    total = 0390    for axis in range(occ.ndim):391        total += int((occ & _shift_occ(occ, axis)).sum())392    return total393394# CORE GRAPH395396def core_nodes(design_fn: Callable, number: int, level: int) -> int:397    dim = design_fn(number).ndim398    fill = _FILL_FN[(design_fn.__name__, dim)](number, level)399    return fill400401def core_edges(design_fn: Callable, number: int, level: int) -> int:402    base = design_fn(number)403    e1 = _core_edges_grid(base)404    if e1 == 0:405        return 0406    if level == 1:407        return e1408    g = base409    seq = [e1]410    for _ in range(3):411        g = np.kron(g, base).astype(np.uint8)412        seq.append(_core_edges_grid(g))413    det = seq[1] * seq[1] - seq[0] * seq[2]414    M = [[seq[1], -seq[0]], [seq[2], -seq[1]]]415    rhs = [seq[2], seq[3]]416    d = M[0][0] * M[1][1] - M[0][1] * M[1][0]417    s = Fraction(rhs[0] * M[1][1] - M[0][1] * rhs[1], d)418    p = Fraction(M[0][0] * rhs[1] - rhs[0] * M[1][0], d)419    out = list(seq)420    while len(out) < level:421        out.append(int(s * out[-1] - p * out[-2]))422    return int(out[level - 1])423424# TUNNEL GRAPH425426def tunnel_nodes(design_fn: Callable, number: int, level: int) -> int:427    dim = design_fn(number).ndim428    grid = _GRID_FN[dim](number, level)429    fill = _FILL_FN[(design_fn.__name__, dim)](number, level)430    return grid - fill431432def _tunnel_edges_grid(grid: np.ndarray) -> int:433    inverted = 1 - (grid != 0).astype(np.uint8)434    return _core_edges_grid(inverted)435436def tunnel_edges(design_fn: Callable, number: int, level: int) -> int:437    grid = design_fn(number)438    base = grid439    for _ in range(1, level):440        grid = np.kron(grid, base).astype(np.uint8)441    return _tunnel_edges_grid(grid)442443def _slice_start(number: int, level: int) -> int:444    return 0445446def _slice_triangle(x: int, y: int, start: int):447    north = (x + y + start) % 2 == 0448    if north:449        return [(x, 2 * y + 2), (x + 1, 2 * y), (x + 2, 2 * y + 2)]450    return [(x, 2 * y), (x + 1, 2 * y + 2), (x + 2, 2 * y)]451452def _slice_edges_of(corners):453    a, b, c = corners454    return [tuple(sorted((a, b))), tuple(sorted((b, c))), tuple(sorted((a, c)))]455456def _slice_cut_types(design_fn: Callable, number: int, level: int):457    from mrlypy.six.geometry import cut458    from mrlypy.three.models import Cell3d459    grid = design_fn(number)460    base = grid461    for _ in range(1, level):462        grid = np.kron(grid, base).astype(np.uint8)463    cell6 = cut(Cell3d(types=grid))464    return cell6._cell.types, int(cell6.start)465466def _slice_adjacency_counts(types: np.ndarray, start: int, value: int):467    height, width = types.shape468    cells = [(x, y) for y in range(height) for x in range(width) if int(types[y, x]) == value]469    edge_to_cells: Dict[Tuple, list] = {}470    for (x, y) in cells:471        for e in _slice_edges_of(_slice_triangle(x, y, start)):472            edge_to_cells.setdefault(e, []).append((x, y))473    edges = sum(1 for shared in edge_to_cells.values() if len(shared) == 2)474    index = {c: i for i, c in enumerate(cells)}475    adj = {i: [] for i in range(len(cells))}476    for shared in edge_to_cells.values():477        if len(shared) == 2:478            a, b = shared479            adj[index[a]].append(index[b])480            adj[index[b]].append(index[a])481    seen = [False] * len(cells)482    comps = 0483    for s in range(len(cells)):484        if seen[s]:485            continue486        comps += 1487        stack = [s]488        seen[s] = True489        while stack:490            u = stack.pop()491            for v in adj[u]:492                if not seen[v]:493                    seen[v] = True494                    stack.append(v)495    return edges, comps496497def slice_core_nodes(design_fn: Callable, number: int, level: int) -> int:498    types, start = _slice_cut_types(design_fn, number, level)499    return int((types == 1).sum())500501def slice_tunnel_nodes(design_fn: Callable, number: int, level: int) -> int:502    types, start = _slice_cut_types(design_fn, number, level)503    return int((types == 0).sum())504505def slice_core_edges(design_fn: Callable, number: int, level: int) -> int:506    types, start = _slice_cut_types(design_fn, number, level)507    edges, _ = _slice_adjacency_counts(types, start, 1)508    return edges509510def slice_tunnel_edges(design_fn: Callable, number: int, level: int) -> int:511    types, start = _slice_cut_types(design_fn, number, level)512    edges, _ = _slice_adjacency_counts(types, start, 0)513    return edges514515def slice_core_components(design_fn: Callable, number: int, level: int) -> int:516    types, start = _slice_cut_types(design_fn, number, level)517    _, comps = _slice_adjacency_counts(types, start, 1)518    return comps519520def solid_slice_core_edges(number: int) -> int:521    if number % 2 == 0:522        raise MrlyError("solid slice closed form is defined for odd number = 2k-1.")523    k = (number + 1) // 2524    return 36 * k ** 2 - 42 * k + 12525526def solid_slice_core_nodes(number: int) -> int:527    if number % 2 == 0:528        raise MrlyError("solid slice closed form is defined for odd number = 2k-1.")529    k = (number + 1) // 2530    return 24 * k ** 2 - 24 * k + 6