from PIL import Image
import numpy as np

# Load the image in grayscale
img = np.array(Image.open("output2g.png").convert("L"))

if img is None:
    print("Failed to load image.")
    exit(1)

height, width = img.shape
print(f"Image dimensions: {height} rows x {width} columns")

# Convert to binary matrix (Black (<128) = 1, White (>=128) = 0)
matrix = (img < 128).astype(np.uint8)

# Check sparsity of columns and rows
row_weights = np.sum(matrix, axis=1)
col_weights = np.sum(matrix, axis=0)

print(f"Row weights: min={np.min(row_weights)}, max={np.max(row_weights)}, avg={np.mean(row_weights):.2f}")
print(f"Col weights: min={np.min(col_weights)}, max={np.max(col_weights)}, avg={np.mean(col_weights):.2f}")

# The matrix is 64x8 blocks of 256x256. 
# Let's verify circulant structure by checking the first row of block 0,0 and see if rows 1..255 are just cyclic shifts!
block_h, block_w = 256, 256
num_blocks_y = height // block_h
num_blocks_x = width // block_w

print(f"Analyzing {num_blocks_y}x{num_blocks_x} blocks of size {block_h}x{block_w}")

is_perfectly_qc = True
block_weights = np.zeros((num_blocks_y, num_blocks_x), dtype=int)

for by in range(num_blocks_y):
    for bx in range(num_blocks_x):
        # Extract block
        block = matrix[by*block_h:(by+1)*block_h, bx*block_w:(bx+1)*block_w]
        
        # Check first row
        first_row = block[0, :]
        weight = np.sum(first_row)
        block_weights[by, bx] = weight
        
        # Verify cyclic property
        for r in range(1, block_h):
            expected_row = np.roll(first_row, r)
            if not np.array_equal(block[r, :], expected_row):
                is_perfectly_qc = False

print(f"Is perfectly Quasi-Cyclic? {is_perfectly_qc}")

print("\nBlock Row-Weight Map (64 rows of 8 blocks):")
for by in range(min(16, num_blocks_y)):  # Print first 16
    print([block_weights[by, bx] for bx in range(num_blocks_x)])
if num_blocks_y > 16:
    print("...")

# Let's check the rank of the matrix over GF(2)!
# We can't do 16384x2048 easily in native numpy, but we can do it via rref or just evaluating over GF(2)
# Actually, since it's 16384x2048, rank is at most 2048. 
# Wait, let's look for linear dependencies or sparse rows directly.

# Are there any suspiciously low-weight blocks?
min_block_weight = np.min(block_weights)
print(f"\nMinimum block weight observed: {min_block_weight}")
if min_block_weight < 64:
    coords = np.where(block_weights == min_block_weight)
    print(f"Low weight blocks located at: {list(zip(coords[0], coords[1]))}")

# Let's analyze the polynomials in the first block row to see if they share a common sparse factor.
