[1]:
import time
import pickle
import math
import numpy as np
from scipy.linalg import solve, LinAlgError
from itertools import combinations
from multiprocess import Pool, cpu_count
from functools import lru_cache
from tqdm import tqdm
from chipsplitting.utils import gauss, get_array_index, to_coordinate

Finding fundamental models

We first reduced our search space by excluding support candidates in \(V_d\), see the support_candidates.ipynb jupyter notebook.

Next, we check a support candidate \(S = \{ (\nu_i, \mu_i) \}_{i=0}^n\) is the support of a fundamental model by solving for the scaling coefficients \(c_i\) that satisfy the constraint \(\sum_{i=0}^n c_i t^{\nu_i} (1-t)^{\mu_i} = 1\). If the scalings are unique and non-negative, we have found a fundamental model.

The following algorithm is a parallelized version of [BM25_MATHREPO], Step 4.

[2]:
@lru_cache(maxsize=None)
def expand_expression(expression):
    if expression.startswith('t^'):
        power = int(expression[2:])
        return tuple([1 if k == power else 0 for k in range(power + 1)])
    elif expression.startswith('(1-t)^'):
        power = int(expression[len('(1-t)^'):])
        return tuple([math.comb(power,k) * (-1)**(k) for k in range(power + 1)])
    else:
        raise ValueError("Invalid expression")

@lru_cache(maxsize=None)
def multiply_expression(left_exp: list[int], right_exp: list[int]):
    left_power = len(left_exp) - 1
    right_power = len(right_exp) - 1
    degree = left_power + right_power
    t_coefficients = [0] * (degree + 1)
    for a in range(left_power + 1):
        for b in range(right_power + 1):
            t_coefficients[a + b] += left_exp[a] * right_exp[b]
    return tuple(t_coefficients)

def pad_with_zero(array, desired_length: int):
    if len(array) < desired_length:
        return list(array) + [0] * (desired_length - len(array))
    return array

def create_matrix_from_support(degree: int, support: list[int]):
    support = [to_coordinate(index) for index in support]
    A = np.array([pad_with_zero(multiply_expression(expand_expression(f"t^{col}"), expand_expression(f"(1-t)^{row}")), degree + 1) for col, row in support])
    return A.T

def is_fundamental(support_as_index, degree):
    A = create_matrix_from_support(degree, support_as_index)
    b = np.array([1] + [0] * degree)
    sol, norm, rank, singular_values = np.linalg.lstsq(A, b, rcond=None)

    if rank < len(support_as_index):
        return None

    if len(norm) == 0:
        return sol
    if np.isclose(norm[0], 0):
        if norm[0] <1e-20:
            return sol
        else:
            print(norm)
            assert False

    return None

def apply_symmetry(config):
    def swap(tuple):
        return (tuple[1], tuple[0])
    def to_array_index(tuple):
        return get_array_index(tuple[0], tuple[1])
    return tuple(sorted([to_array_index(swap(to_coordinate(index))) for index in config]))

def is_positive(solution):
    for x in solution:
        if np.isclose(x,0):
            return False
        if x < 0:
            return False
    return True

def is_parametric(solution):
    for x in solution:
        if not x.is_number:
            return True
    return False

def is_asymmetric(support):
    return list(support) != list(apply_symmetry(support))


def process_single_support(args):
    support, n, d = args
    fundamental = set()

    if len(support) <= n:
        num_selections = n + 1 - len(support)
        domain = [x for x in range(gauss(d + 1)) if x not in support]
        S = [tuple(sorted(support + c)) for c in combinations(domain, num_selections)]
    elif len(support) == n + 1:
        S = [support]
    else:
        return fundamental  # Shouldn't happen

    for s in S:
        solution = is_fundamental(s, d)
        if solution is not None and is_positive(solution):
            fundamental.add(tuple(s))
            if is_asymmetric(s):
                fundamental.add(tuple(sorted(apply_symmetry(s))))

    return fundamental

def find_fundamental_models_parallel(n, d, supports_to_check):
    start = time.time()
    num_workers = cpu_count()
    print(f"Using {num_workers} parallel processes")

    args = [(tuple(support), n, d) for support in supports_to_check]

    fundamental_models = set()
    with Pool(processes=num_workers) as pool:
        for result in tqdm(pool.imap_unordered(process_single_support, args), total=len(args)):
            fundamental_models.update(result)

    end = time.time()

    return fundamental_models

Fundamental models in \(\Delta_7\)

[3]:
n, d = 7, 7

with open(f'data/n{n:02}_d{d:02}_minimal.pkl', 'rb') as f:
    data = [list(candidate) for candidate in pickle.load(f)]
    iter_start_time = time.perf_counter()
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    iter_elapsed_time = time.perf_counter() - iter_start_time
    print(f"n={n:<2}, d={d:<2} | Fundamental models found: {len(fundamental_models):<7} | Time: {iter_elapsed_time:>5.2f}s")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
Using 12 parallel processes
100%|█████████████████████████████████████████████████████████████████| 15779/15779 [00:16<00:00, 967.19it/s]
n=7 , d=7  | Fundamental models found: 83906   | Time: 16.43s

[16]:
n, d = 7, 8

with open(f'data/n{n:02}_d{d:02}_minimal.pkl', 'rb') as f:
    data = [list(candidate) for candidate in pickle.load(f)]
    iter_start_time = time.perf_counter()
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    iter_elapsed_time = time.perf_counter() - iter_start_time
    print(f"n={n:<2}, d={d:<2} | Fundamental models found: {len(fundamental_models):<7} | Time: {iter_elapsed_time:>5.2f}s")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
Using 12 parallel processes
100%|████████████████████████████████████████████████████████████████| 64000/64000 [00:48<00:00, 1316.09it/s]
n=7 , d=8  | Fundamental models found: 23285   | Time: 48.74s

[4]:
n, d = 7, 9

with open(f'data/n{n:02}_d{d:02}.pkl', 'rb') as f:
    data = [list(candidate) for candidate in pickle.load(f)]
    iter_start_time = time.perf_counter()
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    iter_elapsed_time = time.perf_counter() - iter_start_time
    print(f"n={n:<2}, d={d:<2} | Fundamental models found: {len(fundamental_models):<7} | Time: {iter_elapsed_time:>5.2f}s")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
Using 12 parallel processes
100%|████████████████████████████████████████████████████████████| 1927201/1927201 [04:48<00:00, 6690.52it/s]
n=7 , d=9  | Fundamental models found: 6445    | Time: 288.47s

[5]:
n, d = 7, 10
with open(f'data/n{n:02}_d{d:02}.pkl', 'rb') as f:
    data = [list(candidate) for candidate in pickle.load(f)]
    iter_start_time = time.perf_counter()
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    iter_elapsed_time = time.perf_counter() - iter_start_time
    print(f"n={n:<2}, d={d:<2} | Fundamental models found: {len(fundamental_models):<7} | Time: {iter_elapsed_time:>5.2f}s")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
Using 12 parallel processes
100%|████████████████████████████████████████████████████████████| 6508380/6508380 [14:32<00:00, 7459.34it/s]
n=7 , d=10 | Fundamental models found: 1442    | Time: 873.87s
[3]:
n, d = 7, 11
with open(f'data/n{n:02}_d{d:02}.pkl', 'rb') as f:
    data = [list(candidate) for candidate in pickle.load(f)]
    iter_start_time = time.perf_counter()
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    iter_elapsed_time = time.perf_counter() - iter_start_time
    print(f"n={n:<2}, d={d:<2} | Fundamental models found: {len(fundamental_models):<7} | Time: {iter_elapsed_time:>5.2f}s")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
Using 12 parallel processes
100%|██████████████████████████████████████████████████████████| 18676991/18676991 [39:50<00:00, 7813.80it/s]
n=7 , d=11 | Fundamental models found: 332     | Time: 2393.49s
[4]:
n, d = 7, 12
with open(f'data/n{n:02}_d{d:02}.pkl', 'rb') as f:
    data = [list(candidate) for candidate in pickle.load(f)]
    iter_start_time = time.perf_counter()
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    iter_elapsed_time = time.perf_counter() - iter_start_time
    print(f"n={n:<2}, d={d:<2} | Fundamental models found: {len(fundamental_models):<7} | Time: {iter_elapsed_time:>5.2f}s")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
Using 12 parallel processes
100%|████████████████████████████████████████████████████████| 47682475/47682475 [1:38:41<00:00, 8051.82it/s]
n=7 , d=12 | Fundamental models found: 56      | Time: 5930.48s
[4]:
n, d = 7, 13
with open(f'data/n{n:02}_d{d:02}.pkl', 'rb') as f:
    data = [list(candidate) for candidate in pickle.load(f)]
    iter_start_time = time.perf_counter()
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    iter_elapsed_time = time.perf_counter() - iter_start_time
    print(f"n={n:<2}, d={d:<2} | Fundamental models found: {len(fundamental_models):<7} | Time: {iter_elapsed_time:>5.2f}s")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
n=7 | d=13 | 107547676 supports to check (parallelized)
Using 12 parallel processes
100%|██████████████████████████████████████████████████████| 107547676/107547676 [3:53:50<00:00, 7665.30it/s]
Solved linear equations for all supports. Elapsed time: 14036.57s

We have found 8 fundamental models.
There exist exactly 8 unique fundamental models ✨
Serialized to fundamental-models/n07_d13.pkl
------------------------------------

\(n = 8\)

[4]:
n, d = 8, 8
with open(f'data/n{n:02}_d{d:02}_minimal.pkl', 'rb') as f:
    data = pickle.load(f)
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    print(f"There exist exactly {len(fundamental_models)} unique fundamental models ✨")
    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
        print(f"Serialized to {f.name}")
        print("------------------------------------")
n=8 | d=8 | 66830 supports to check (parallelized)
Using 12 parallel processes
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 66830/66830 [11:33<00:00, 96.41it/s]
Solved linear equations for all supports. Elapsed time: 693.29s

We have found 1279349 fundamental models.
There exist exactly 1279349 unique fundamental models ✨
Serialized to fundamental-models/n08_d08.pkl
------------------------------------
[4]:
n, d = 8, 9
with open(f'data/n{n:02}_d{d:02}_minimal.pkl', 'rb') as f:
    data = pickle.load(f)
    fundamental_models = find_fundamental_models_parallel(n, d, data)
    print(f"There exist exactly {len(fundamental_models)} unique fundamental models ✨")

    with open(f'fundamental-models/n{n:02d}_d{d:02d}.pkl', 'wb') as f:
        pickle.dump(fundamental_models, f)
        print(f"Serialized to {f.name}")
        print("------------------------------------")
n=8 | d=9 | 462075 supports to check (parallelized)
Using 12 parallel processes
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 462075/462075 [43:08<00:00, 178.52it/s]
Solved linear equations for all supports. Elapsed time: 2588.41s

We have found 285601 fundamental models.
There exist exactly 285601 unique fundamental models ✨
Serialized to fundamental-models/n08_d09.pkl
------------------------------------

Visualizing fundamental models in \(\Delta_7\) of degree \(13\)

We draw all fundamental models in \(\Delta_7\) of degree \(13\).

[32]:
import chipsplitting as cs
n, d = 7, 13
with open(f'fundamental-models/n{n:02}_d{d:02}.pkl', 'rb') as f:
    fundamental_models = pickle.load(f)
    for m in fundamental_models:
        print([to_coordinate(x) for x in m])
        l = cs.LinearForm(np.array([1 if i in m else 0 for i in range(gauss(d+1))]), np.array([0 if i > 0 else 1 for i in range(gauss(d+1))])).to_hyperfield()
        print(l)
        print("----------------------------------------")
[(1, 6), (3, 5), (5, 4), (7, 3), (9, 2), (11, 1), (0, 13), (13, 0)]
   +
   .    .
   .    .    .
   .    .    .    .
   .    .    .    .    .
   .    .    .    .    .    .
   .    .    .    .    .    .    .
   .    +    .    .    .    .    .    .
   .    .    .    +    .    .    .    .    .
   .    .    .    .    .    +    .    .    .    .
   .    .    .    .    .    .    .    +    .    .    .
   .    .    .    .    .    .    .    .    .    +    .    .
   .    .    .    .    .    .    .    .    .    .    .    +    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------
[(6, 1), (5, 3), (4, 5), (3, 7), (2, 9), (1, 11), (0, 13), (13, 0)]
   +
   .    .
   .    +    .
   .    .    .    .
   .    .    +    .    .
   .    .    .    .    .    .
   .    .    .    +    .    .    .
   .    .    .    .    .    .    .    .
   .    .    .    .    +    .    .    .    .
   .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    +    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    +    .    .    .    .    .    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------
[(1, 1), (1, 6), (8, 2), (7, 4), (1, 11), (11, 1), (0, 13), (13, 0)]
   +
   .    .
   .    +    .
   .    .    .    .
   .    .    .    .    .
   .    .    .    .    .    .
   .    .    .    .    .    .    .
   .    +    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    +    .    .
   .    .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    +    .    .    .
   .    +    .    .    .    .    .    .    .    .    .    +    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------
[(1, 1), (6, 1), (2, 8), (4, 7), (1, 11), (11, 1), (0, 13), (13, 0)]
   +
   .    .
   .    +    .
   .    .    .    .
   .    .    .    .    .
   .    .    +    .    .    .
   .    .    .    .    +    .    .
   .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .    .    .
   .    +    .    .    .    .    +    .    .    .    .    +    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------
[(3, 3), (6, 1), (3, 7), (7, 3), (2, 9), (1, 11), (0, 13), (13, 0)]
   +
   .    .
   .    +    .
   .    .    .    .
   .    .    +    .    .
   .    .    .    .    .    .
   .    .    .    +    .    .    .
   .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .
   .    .    .    +    .    .    .    +    .    .    .
   .    .    .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    +    .    .    .    .    .    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------
[(1, 1), (1, 6), (7, 3), (9, 2), (1, 11), (11, 1), (0, 13), (13, 0)]
   +
   .    .
   .    +    .
   .    .    .    .
   .    .    .    .    .
   .    .    .    .    .    .
   .    .    .    .    .    .    .
   .    +    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    +    .    .    .
   .    .    .    .    .    .    .    .    .    +    .    .
   .    +    .    .    .    .    .    .    .    .    .    +    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------
[(1, 1), (6, 1), (3, 7), (2, 9), (1, 11), (11, 1), (0, 13), (13, 0)]
   +
   .    .
   .    +    .
   .    .    .    .
   .    .    +    .    .
   .    .    .    .    .    .
   .    .    .    +    .    .    .
   .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .    .    .
   .    +    .    .    .    .    +    .    .    .    .    +    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------
[(3, 3), (1, 6), (3, 7), (7, 3), (9, 2), (11, 1), (0, 13), (13, 0)]
   +
   .    .
   .    .    .
   .    .    .    .
   .    .    .    .    .
   .    .    .    .    .    .
   .    .    .    +    .    .    .
   .    +    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .
   .    .    .    .    .    .    .    .    .    .
   .    .    .    +    .    .    .    +    .    .    .
   .    .    .    .    .    .    .    .    .    +    .    .
   .    .    .    .    .    .    .    .    .    .    .    +    .
   -    .    .    .    .    .    .    .    .    .    .    .    .    +

----------------------------------------

References

[BM25] Arthur Bik and Orlando Marigliano. Classifying one-dimensional discrete models with maximum likelihood degree one. Advances in Applied Mathematics, 170:102928, 2025.

[BM25_MATHREPO] https://mathrepo.mis.mpg.de/ChipsplittingModels/Section8.html

[ ]: