Order elements in a matrix based on another matrix

Question

I'd like to order an arbitrary number of elements into an arbitrarily shaped matrix (matrix_a) based on values in a binary matrix (element_map) that describes attributes of the elements. matrix_b defines which elements can be adjacent to each other in matrix_a. "Adjacent" in this case includes diagonal. A concrete, but toy example:

import numpy as np

#This is a just a mask of the final matrix, to show the shape.
#The center position (5) is adjacent to all other positions whereas
#position (1) is adjacent to 4, 2 and 5, etc. 
matrix_a = np.array([[1, 4, 7],
                     [2, 5, 8],
                     [3, 6, 9]])

elements = range(0, 10)

#each 'row' corresponds to an element
#each 'column' corresponds to an attribute of the various elements
#here, each element has 5 attributes, which are always 0 or 1
element_map = np.array([[0, 0, 1, 0, 1], #element 1
                        [1, 1, 1, 0, 1], #element 2
                        [1, 0, 0, 1, 0], #element 3
                        [1, 0, 1, 0, 0], #etc.
                        [1, 1, 1, 0, 0],
                        [1, 1, 0, 1, 0],
                        [1, 0, 1, 0, 1],
                        [1, 1, 0, 0, 0],
                        [1, 1, 0, 0, 0]]) #element 9

The desired output is that the nine elements are placed in matrix_a, each only appearing once, based on an adjacency rule for which matrix_b is needed. The rule is that for any element placed in matrix_a, the element's value for attribute x (can be any attribute) must be zero, and the value of all adjacent (adjacent in matrix_a) elements must be 1 for attribute x.

Assume also that we have a function that can accept a matrix and coordinates and return all the adjacent values:

def adj(m, x, y):
    ''' find adjacent values to coordinates x,y in a matrix, m'''
    if x == 0:
        if y == 0:
            r = m[0:2,0:2] # x==0, y==0
        else:
            r = m[0:2,y-1:y+2]  # x==0, y!=0
    elif y == 0: # x != 0, y == 0
        r = m[x-1:x+2,0:2]
    else: #any other position
        r = m[(x-1):(x+2),(y-1):(y+2)]
        
    return [i for i in r.flatten().tolist() if i != m[x,y]]

#example:
q = np.arange(1,97).reshape(12,8).T

array([[ 1,  9, 17, 25, 33, 41, 49, 57, 65, 73, 81, 89],
       [ 2, 10, 18, 26, 34, 42, 50, 58, 66, 74, 82, 90],
       [ 3, 11, 19, 27, 35, 43, 51, 59, 67, 75, 83, 91],
       [ 4, 12, 20, 28, 36, 44, 52, 60, 68, 76, 84, 92],
       [ 5, 13, 21, 29, 37, 45, 53, 61, 69, 77, 85, 93],
       [ 6, 14, 22, 30, 38, 46, 54, 62, 70, 78, 86, 94],
       [ 7, 15, 23, 31, 39, 47, 55, 63, 71, 79, 87, 95],
       [ 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96]])

adj(q, 5, 5)
[37, 45, 53, 38, 54, 39, 47, 55]

As a side note, the possible search space for an 8x12 matrix and 96 elements is 10**149 = 96!

Maybe this would be a job for linear/integer programming, but I'm not sure how to formulate the adjacency constraints

Answer 1

This can be done with an ILP (Integer Linear Program). A rough cut is shown below. It seems that you have no weighting declared to consider one result superior to another, so you might also have luck with a CP (Constraint Programming) approach. The ILP below uses the number of successful element fills as an objective. There may be many solutions (or none), and it will just grab the first discovered.

Of note, the example you provided has no solution for a 3x3 target matrix. I used the segment below that to make elements from the identity matrix (flipped 1/0), which--by inspection--is the best at fitting element list possible.

I tweaked your adj() function a bit because you had a mixed schema where you were taking in (x, y) and returning a numbered position.... I just made it all (x, y)

If there is a solution, this tends to find it very quickly for moderate sizes of matrix_a, especially if the solution set is large.

I used cbc solver here, which is nice freeware, but there are other MIP solvers that work fine (glpk, highs, ...)

There is still some "meat on the bone" here to optimize this a bit, mostly by domain trimming the [e, p, a] tuples down to legal assignments based on the a value, meaning any place[e, p, a] where element e has a non-zero for attribute a is not legal and should be removed from the domain for place. Then you might look at pre-computing the possible neighbors for any given element to reduce the size of the adjacency constraint a bit.

Code:

"""
fit elements into matrix based on adjacency rules
"""

import numpy as np
import pyomo.environ as pyo


class Element:
    """a convenience to hold the rows of attribute values"""
    def __init__(self, row):
        self.attributes = tuple(row)

    def attribute(self, idx):
        return self.attributes[idx]

    def __repr__(self):
        return str(self.attributes)


class Position:
    """a convenience for (x, y) positions that must have equality & hash defined for consistency"""
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def __repr__(self):
        return f'({self.x}, {self.y})'

    def __hash__(self):
        return hash((self.x, self.y))

    def __eq__(self, other):
        if isinstance(other, Position):
            return (self.x, self.y) == (other.x, other.y)
        return False

# each 'row' corresponds to an element
# each 'column' corresponds to an attribute of the various elements
# here, each element has 5 attributes, which are always 0 or 1
element_map = np.array([[0, 0, 1, 0, 1],  # element 1
                        [1, 1, 1, 0, 1],  # element 2
                        [1, 0, 0, 1, 0],  # element 3
                        [1, 0, 1, 0, 0],  # etc.
                        [1, 1, 1, 0, 0],
                        [1, 1, 0, 1, 0],
                        [1, 0, 1, 0, 1],
                        [1, 1, 0, 0, 0],
                        [1, 1, 0, 0, 0]])  # element 9

# Alternate element map...
d = 5
element_map = np.ones((d, d), dtype=int) - np.eye(d, dtype=int)
print(element_map)

matrix_a_rows = 3
matrix_a_cols = 3
matrix_a = np.zeros((matrix_a_rows, matrix_a_cols))


def adj_xy(mat, p: Position):
    x, y = p.x, p.y
    res = []
    rows = len(mat) - 1
    cols = len(mat[0]) - 1
    for i in range(x - 1, x + 2):
        for j in range(y - 1, y + 2):
            if all((0 <= i <= rows, 0 <= j <= cols, (i, j) != (x, y))):
                res.append(Position(i, j))
    return res


# SET UP ILP
m = pyo.ConcreteModel('matrix_fitter')

# SETS

m.E = pyo.Set(initialize=[Element(row) for row in element_map], doc='elements')
m.P = pyo.Set(initialize=[Position(x, y) for x in range(len(matrix_a)) for y in range(len(matrix_a[0]))],
              doc='positions')
m.A = pyo.Set(initialize=list(range(len(element_map[0]))), doc='attribute')

# VARS
# place element e in position p based on attribute a being 0...
m.place = pyo.Var(m.E, m.P, m.A, domain=pyo.Binary, doc='place')

# OBJ:  Place as many as possible....  ideally --> len(matrix_a)
m.obj = pyo.Objective(expr=pyo.sum_product(m.place), sense=pyo.maximize)


# CONSTRAINTS
@m.Constraint(m.E, m.P, m.A)
def zero_attribute(m, e: Element, p, a):
    # can only place if attribute a is zero
    return m.place[e, p, a] * e.attribute(a) == 0


@m.Constraint(m.P)
def only_one(m, p):
    # can only place at most one element at each location
    return sum(m.place[e, p, a] for e in m.E for a in m.A) <= 1


@m.Constraint(m.E, m.P, m.A)
def adjacency(m, e: Element, p: Position, a):
    # neighbors must "add up" in the selected attribute (a), if placed in neighbor positions
    # we can "add up" the values of the selected attribute (a) for all adjacent positions that were selected
    neighbor_positions = adj_xy(matrix_a, p)
    return (sum(m.place[ee, pp, aa] * ee.attribute(a) for ee in m.E for pp in neighbor_positions for aa in m.A) >=
            len(neighbor_positions) * m.place[e, p, a])


solver = pyo.SolverFactory('cbc')
results = solver.solve(m, tee=True)
print(results)
if results.solver.termination_condition == pyo.TerminationCondition.optimal:
    for idx in m.place.index_set():
        if m.place[idx].value == 1:
            print(idx)
    if pyo.value(m.obj) == matrix_a_rows * matrix_a_cols:
        # all positions were filled
        print('success!')
    else:
        print(f'the max number of elements that can be placed is {pyo.value(m.obj)} / {matrix_a_rows * matrix_a_cols}')

else:
    print('Problem with model...see results')

Output (truncated):

((0, 1, 1, 1, 1), (0, 1), 0)
((0, 1, 1, 1, 1), (2, 1), 0)
((1, 1, 0, 1, 1), (1, 0), 2)
((1, 1, 0, 1, 1), (1, 2), 2)
((1, 1, 1, 0, 1), (1, 1), 3)
((1, 1, 1, 1, 0), (0, 0), 4)
((1, 1, 1, 1, 0), (0, 2), 4)
((1, 1, 1, 1, 0), (2, 0), 4)
((1, 1, 1, 1, 0), (2, 2), 4)
success!