Formulating a constraint with pyomo

Question

I'm trying to array genes (elements) in a matrix (matrix_a) based on the restriction sites they have (attributes). Restriction sites are recognized by enzymes that cut the DNA at a specific site. This can be represented as a sparse matrix with rows as genes, and columns as restriction enzymes, and values being zero or one, meaning the enzyme cuts the gene or it doesn't cut the gene. Given that information, I want to use ILP to array the genes such that the number of enzymes needed for each gene is minimized.

Background: Given a 3X3 matrix, I'm trying to arrange genes (elements) in the matrix (matrix_a) based on an adjacency rule. For any given position in the matrix, a restriction enzyme (attribute) should be used that cuts genes in all adjacent positions in the matrix, but not the gene at that position. I got help to produce an adjacency constraint that works. I just wanted to provide an overall context for the problem.

The problem is that sometimes it'll be necessary to add more than one restriction enzyme to a gene in order to cut all of the neighboring genes, and I'm stuck trying to formulate this as a constraint. My attempt at doing this immediately follows the CONSTRAINTS section in the code below; the code is commented-out. Ideally, the constraints would be:

1. Only one gene (element) can be assigned to a position in the matrix, but within that one position, the gene (element) can appear multiple times if multiple restriction enzymes (attributes) are used
2. Once a gene (element) is assigned to a position in the matrix, it must not appear in any other positions

One way of doing this is: when I iterate over all elements and all positions, I need to collapse cases where one element appears at a given position more than once into a single value rather than the taking the sum. I haven't been able to distinguish, via constraints, cases where a gene (element) appears in a single position more than once and cases where a gene (element) appears in multiple positions.

Note I'm not sure if the element_map in my code below actually has a solution given my desired constraints

import pyomo.environ as pyo
import numpy as np

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


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([[1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0],
                        [1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1],
                        [1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1],
                        [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
                        [0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1],
                        [1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
                        [1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1],
                        [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1],
                        [1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1]])

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')



#only consider valid element positionings
valid_vals = []
for element in m.E:
    for pos in m.P:
        for i in m.A:
            if element.attribute(i) == 0:
                valid_vals.append((element,pos,i))

m.VS = pyo.Set(initialize=set(valid_vals))

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


# OBJ:  
m.obj = pyo.Objective(expr=pyo.sum_product(m.place), sense=pyo.minimize)

# CONSTRAINTS

#each element can only appear in one position, but can appear multiple times inthat position with different attributes
#sum across all positions (not within one position) must be == 1
#collapse cases where element is in one position twice to one
#something where 0==0, 1==1, 2==1
#m.one_element_per_across_all_positions = pyo.ConstraintList()
#for e in m.E:
#    for p in m.P:
#        s = 0
#        for a in m.A:
#            if e.attribute(a) == 1: continue
#            
#            s += m.place[e,p,a]
#       
#        m.one_element_per_across_all_positions.add(s <= 2)

#each place must have 1 or more elements
m.single = pyo.ConstraintList()
for p in m.P:
    s = 0
    for e in m.E:
        for a in m.A:
            if e.attribute(a) == 1: continue
            s += m.place[e,p,a]
    m.single.add(s >= 1)

#each element/attribute combo appears only once
m.eacombo = pyo.ConstraintList()
for element in m.E:
    for attribute in m.A:
        if element.attribute(attribute) == 1: continue
        s = 0
        for p in m.P:
            s += m.place[element,p,attribute]
        m.eacombo.add(s <= 1)

#adjacency constraint
m.attr_constraint = pyo.ConstraintList()
for p in m.P:
    for e in m.E:
        s = 0
        for a in m.A:
            if e.attribute(a) == 1: continue
            s += m.place[e,p,a]
        m.attr_constraint.add(s <= 2)

m.adjacency_rule = pyo.ConstraintList()
for place in m.P:
    neighbor_positions = adj_xy(matrix_a, place)
    for element in m.E:
        for attribute in m.A:
            if element.attribute(attribute) == 1 : continue
            s = 0
            for ee in m.E:
                if ee.attribute(attribute) == 0: continue
                if ee == element: continue
               
                for pp in neighbor_positions:
                    for aa in m.A:
                        if ee.attribute(aa) == 1: continue
                        s += m.place[ee,pp,aa] * ee.attribute(attribute)

            m.adjacency_rule.add(s >= len(neighbor_positions)*m.place[element,place,attribute])


solver = pyo.SolverFactory('cbc')

results = solver.solve(m, tee=True)
print(results, 'results')
if results.solver.termination_condition == pyo.TerminationCondition.optimal:
    for idx in m.place.index_set():
        if m.place[idx].value == 1:
            print(idx, 'idx')
    if pyo.value(m.obj) == matrix_a_rows * matrix_a_cols:
        # all positions were filled
        print('success!')
    else:
        print(f'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')

The result I get is:

((1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1), (2, 2), 1) idx
((1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1), (1, 1), 3) idx
((1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1), (2, 0), 1) idx
((0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1), (1, 0), 0) idx
((1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1), (0, 0), 1) idx
((0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1), (1, 2), 11) idx
((1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0), (2, 1), 5) idx
((1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1), (0, 2), 1) idx
((1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0), (0, 1), 12) idx

Notice that the element in position (2,1) is the same as the element in (0,1). This would be equivalent to having the same gene in two different positions in the matrix, which is what I don't want. I need a constraint that a) prevents the same element from appearing in two different positions and also b) allows for the same element to appear within one position multiple times, with different attributes