# Data structure for finding pairs of matrices under constraint

PROBLEM DESCRIPTION

Given three positive integers k, n1, n2, I have a list of matrices

matrices = Tuples[Range[0,k], {n1,n2}]


Since Length[matrices] = $(k+1)^{n1 \cdot n2}$ grows pretty quickly, we can consider k = 2, n1 = 2, n2 = 3 here.

We can compute the row- and column sums of the matrices:

rowsums = Total[matrices, {2}];
colsums = Total[matrices, {3}];


Now we are given two more non-negative integers p1, p2, e.g. p1 = 0, p2 = 1. I want to find all pairs of matrices {m1, m2} from matrices that satisfy the following constraints:

1. The row sums of m1 - m2 are sorted; Sort[Total[m1-m2, {2}]] == Total[m1 - m2, {2}].
2. The smallest row sum is greater or equal to -p1.
3. The column sums of m1 - m2 are sorted.
4. The smallest column sum is greater or equal to -p2.

Notice in particular that m1 = m2 satisfies the constraints.

Looping through all possible pairs and checking whether or not they satisfy the constrains quickly becomes very time consuming. My idea to solve this without looping through all the pairs, is to put matrices into some type of data structure that would let me find valid m1 for a given m2.

QUESTION: How can we make a data structure suited for this in Mathematica?

My attempt is the following:

tree[_][__] = {};
Do[

Do[
treeval = Take[colsums[[ind]], i];
tree[i][treeval] = Union[tree[i][treeval], {colsums[[ind, i + 1]]}]
, {i, n1 - 1}
];
treeval = colsums[[ind]];
tree[n1][treeval] = Union[tree[n1][treeval], {rowsums[[ind, 1]]}];
Do[
treeval = Take[rowsums[[ind]], j];
tree[n1 + j][colsums[[ind]], treeval] =
Union[tree[n1 + j][colsums[[ind]], treeval], {rowsums[[ind, j + 1]]}]
, {j, n2 - 1}
];

, {ind, Length[matrices]}
];


This works in the following way: matrices are sorted first on their colsums, then on their rowsums. tree[ind][colsumSubset] returns a list of possible values of the ind + 1th colsum, given that the ind first colsums is given by colsumSubset. If the returned list is empty, that means no such matrix exists in matrices. If all colsums have been specified, the next form is tree[n1 + ind - 1][colsums, rowsumSubset]. But I'm note sure this type of construction is anywhere close to ideal for this problem...

Any insights much appreciated!