slicing the hypercube, a set cover problem

slicing_the_hypercube.py

import numpy as np
import cvxpy as cp
import itertools


def collect_vertices(n):
    if n <= 0:
        return np.array([])

    # Generate all possible combinations of +1 and -1 using meshgrid
    meshgrid_values = np.meshgrid(*([-1, 1] for _ in range(n)), indexing='ij')
    
    # Stack the values along the last axis to get the final array
    sign_vectors = np.stack(meshgrid_values, axis=-1).reshape(-1, n)

    return sign_vectors



def collect_vertices_slow(n):
    # Generate all possible combinations of +1 and -1
    sign_vectors = np.array(list(itertools.product([-1, 1], repeat=n)))

    return sign_vectors


def collect_lines(n):
    vertices = collect_vertices(n)
    print(len(vertices))
    line_directions = []
    lines = []
    for i in range(n):
        vs = vertices[vertices[:, i] == -1]
        vs[:, i] = 0
        line_directions.append(i * np.ones(len(vs)))
        lines.append(vs)
    lines = np.array(lines).reshape(-1, n)
    line_directions = np.array(line_directions).flatten()
    print(lines.shape, line_directions.shape)
    return lines, line_directions


def collect_edges(n):
    vertices = collect_vertices(n)
    print(len(vertices))
    v1 = []
    v2 = []
    for i in range(n):
        v1i = vertices[vertices[:, i] == -1]
        v2i = v1i.copy()
        v2i[:, i] = 1
        v1.append(v1i)
        v2.append(v2i)
    v1 = np.array(v1).reshape(-1, n)
    v2 = np.array(v2).reshape(-1, n)
    return v1, v2


def sample_slices(edges, N):
    v1, v2 = edges
    dirs = np.random.normal(size=(N, n))
    biases = np.random.normal(size=(N, ))
    prod1 = dirs @ v1.T + biases[:, None]
    prod2 = dirs @ v2.T + biases[:, None]
    intersects = prod1 * prod2 < 0
    intersects, counts = np.unique(intersects, axis=0, return_counts=True)
    return intersects


n = 6

# lines, line_directions = collect_lines(n)
edges = collect_edges(n)

intersects = sample_slices(edges, N=10000)

plane_count, edge_count = intersects.shape
print("found", plane_count, "distinct slices, slicing", edge_count, "edges")
print(np.histogram(intersects.sum(axis=1)))
largest_slice = intersects.sum(axis=1).max()
print("largest slice size", largest_slice)
print("number of distinct largest slices", (intersects.sum(axis=1) == largest_slice).sum())

# x is a binary variable array where each element represents whether a row is chosen (1) or not (0)
x = cp.Variable(plane_count, boolean=True)

# Objective: Minimize the number of rows chosen
objective = cp.Minimize(cp.sum(x))

# Constraints: The sum of chosen rows should cover each column at most once
constraints = [x @ intersects >= np.ones(edge_count)]

# Define and solve the problem
problem = cp.Problem(objective, constraints)
problem.solve(verbose=True)

# Get the result
chosen_rows = intersects[x.value == 1, :].astype(int)

print("Chosen rows")
print(chosen_rows)