Refining subset relations

People who do combinatorics (like me) are often faced with the following problem:

For a list of combinatorial objects (vectors, permutations graphs,etc.), we know what multi-set of values it generates under some unknown map we wish to explore.

For example: perhaps we know that the four vectors $\{123,213,231,312\}$ generates the list $\{1,2,0,2\}$, but we don't know which permutation give rise to which value. However, if we have several such subset-relations, one can sometimes deduce exactly what value an object gives. Perhaps we also know $\{123, 321 \}$ gives $\{3,0\}$. Since $123$ are common elements in both, and $0$ is the only common statistic, we know that $123$ is responsible for the $0$.

Sometimes, it is not possible to deduce such relations uniquely, but if we have the multiset correspondences $S_1 \to M_1$ and $S_2 \to M_2$, and $S_1 \subseteq S_2$, then we have the refined correspondence $S_1 \cap S_2 \to M_1 \setminus M_2$ as multisets.

Question: Given a list of such multi-set correspondences, produce a new list of maximally refined multiset relations. We also wish to make this minimal, in the sense that if $S_1 \to M_1$ and $S_2 \to M_2$ is in the list, then $S_1 \cup S_2 \to M_1 \cup M_2$ is not in the output, since this is an implicit weaker relation.

EDIT: Below is my implementation so far, but I suspect it can be improved a lot.

The basic idea is to store the multiset of associated values as a sparse vector, since this makes some operation easier to perform.

MultiSetToSparseArray[ms_List, dim_] :=
SparseArray[#1 -> #2 & @@@ Tally[ms], {dim}];
SparseArrayToMultiSet[spA_] :=
Join @@ Table[
ConstantArray[r[[1, 1]], r[]], {r, Most[ArrayRules[spA]]}];

RefineSubsetsRelations[subsetRelations_List] :=
Module[{i, j, s1, s2, m1, m2, int, com,
newRel, foundNew = True, relations = subsetRelations, sa, dim},

(* Maximal possible statistic. *)
dim = Max[Last /@ relations];

(* Sort input subsets. *)

relations = {Sort@#1, MultiSetToSparseArray[#2, dim]} & @@@
relations;

While[ foundNew == True,

foundNew = False;
(* Sort by size. *)

relations = Select[relations, Length[#1] > 0 &];
relations = SortBy[Union@relations, Length[#[]] &];

Do[

If[relations[[i]] === relations[[j]] ,
Continue[]
];

{s1, m1} = relations[[i]];
{s2, m2} = relations[[j]];
int = Intersection[s1, s2];

If[ Length[s1] == 0 || Length[s2] == 0,
Continue[]
];

If[ Length[int] == Length[s1],

If[Min[m2 - m1] < 0,
Complement[s2, s1], ArrayRules[m2 - m1]];
Abort[];
];
(* Replace the more general rule with a specific one. *)

relations[[j]] = {Complement[s2, s1], m2 - m1};

foundNew = True;
Continue[]
];

(* Only one common value,
this means that the intersection only takes this value. *)

com = Intersection[SparseArrayToMultiSet[m1],
SparseArrayToMultiSet[m2]];
If[ Length[com] == 1 && Length[int] > 0,
com = First[com];
sa = SparseArray[com -> Length[int], dim];
newRel = {int, sa};
relations[[i]] = {Complement[s1, int], m1 - sa};
relations[[j]] = {Complement[s2, int], m2 - sa};

foundNew = True;
Continue[]
];

, {i, Length@relations}, {j, i + 1, Length@relations}];
];

{#1, SparseArrayToMultiSet[#2]} & @@@ relations
];


EDIT 2: As for motivation, and lots of examples, see the web site findstat.org This is used to look in a database for combinatorial maps that appear in literature. My example above is borrowed for the number of inversions of a permutation. However, in order to use findstat.org, we need to know which value to assign to each combinatorial object. However, in many research problems, we can only obtain data as above. A famous example is the Kostka-Foulkes polynoials, which basically encodes exactly such a mapping from a subset of combinatorial objects (semi-standard Young tableaux) to a multiset of non-negative integers. By some trickery and extra assumptions, one can get several such intersecting subsets. Once one have some singletons, it is possible to start making conjectures of what the map actually is.

EXAMPLE INPUT Here are 10 equations, where each subset is of size 3, (objects are permutations, or plain vectors if you like), and the multi-set of values is given as the second entry in the pair {setOfPermutations,multiSetOfvalues}.

{
{{{{3, 1, 2}, {1, 3, 2}, {2, 1, 3}}, {1, 1, 2}}},
{{{{3, 2, 1}, {3, 1, 2}, {1, 3, 2}}, {1, 2, 3}}},
{{{{1, 3, 2}, {3, 1, 2}, {2, 3, 1}}, {1, 2, 2}}},
{{{{3, 2, 1}, {2, 3, 1}, {1, 3, 2}}, {1, 2, 3}}},
{{{{1, 3, 2}, {3, 2, 1}, {1, 2, 3}}, {0, 1, 3}}},
{{{{3, 2, 1}, {1, 3, 2}, {2, 1, 3}}, {1, 1, 3}}},
{{{{1, 3, 2}, {3, 2, 1}, {3, 1, 2}}, {1, 2, 3}}},
{{{{2, 3, 1}, {2, 1, 3}, {1, 2, 3}}, {0, 1, 2}}},
{{{{3, 2, 1}, {2, 3, 1}, {3, 1, 2}}, {2, 2, 3}}},
{{{{1, 3, 2}, {2, 3, 1}, {3, 1, 2}}, {1, 2, 2}}}
}

• Are you using "cycle" notation wherein 123 and 231 and 312 are equivalent permutations? If not it might be helpful to have the convention clarified. – Daniel Lichtblau Nov 12 '15 at 20:35
• We are voting to close because we don't understand what you want from us. This is not a code review and improvement service. We will offer advice on code improvement when a question makes it clear what parts of the code are unsatisfactory and why and the question is properly tagged. – m_goldberg Nov 13 '15 at 2:35
• I'm not happy with the wording but I did vote to reopen. The reasons are that (1) It has, after rewording and clarifications, met some level of minimal clarity (maybe emphasis on minimal). (2) It is very much a nontrivial problem. – Daniel Lichtblau Nov 13 '15 at 16:16
• Please add a not-too-large test example. – Daniel Lichtblau Nov 16 '15 at 19:47
• Although I do not really understand your notation, a look at the code shows that it is structurally very similar to what I proposed. A few refinements which I should add are special case handling of singletons and sorting by size with smaller sets considered earlier (I think that would be better for speed but I'm not certain). – Daniel Lichtblau Nov 17 '15 at 0:29

I have some sort of a brute force solution with combinatorial optimization for the test problem. Note that the solution does not use all permutations, only the ones for which the permutations-to-values relationships are known. For the general problem the constraints have to be modified to accomodate the union and complement of the permutations. (I need to think about the problem further...)

The general idea is to use variables x[i,j] that have values 0 or 1. The value of x[i,j] is 1 if to the i-th permutation is assigned the j-th value. We form constraints based on the given equations, and for a fixed i0 we allow only one x[i0,j] to be 1.

I am using Maximize \ Minimize below, but the optimization function and the constriants can be easily recast into a vector and a matrix form respectively. (I.e. amenable for LinearProgramming. See my answer to "How to fill a grid make its total be largest".)

First we map all known permutations into indecies and make the corresponding replacement rules.

perms = Union[Flatten[permsToMultiSetVals[[All, 1]], 1]];
indexToPermRules = Reverse /@ permToIndexRules;
permToIndexRules We do similar mapping for the values.

msVals = Union[Flatten[permsToMultiSetVals[[All, 2]]]];
indexToMsValRules = Reverse /@ msValToIndexRules;
msValToIndexRules The following function makes constraints from one line of the given permutations-to-values relationships.

Clear[MakeConstraints]
MakeConstraints[{perms_, vals_}, varName_Symbol] :=
MakeConstraints[perms, vals, varName];
MakeConstraints[perms_, vals_, varName_Symbol] :=
Block[{t = Tally[vals], pInds, vars},
t[[All, 1]] = t[[All, 1]] /. msValToIndexRules;
pInds = perms /. permToIndexRules;
vars = Outer[varName[#2, #1] &, t[[All, 1]], pInds];
MapThread[Total[#2] == #1 &, {t[[All, 2]], vars}]
];


Example of constraints made for one of the known permutations-to-values relationships:

permsToMultiSetVals[]
MakeConstraints[permsToMultiSetVals[], x] Equations made from all of the known relationships:

eqs = Flatten[MakeConstraints[#, x] & /@ permsToMultiSetVals];
ColumnForm[eqs] All created variables:

vars = Union[Cases[eqs, x[___], \[Infinity]]] Unique mapping per permutation constraints:

uniqueConstraints =
Map[Total[Cases[vars, x[#, _]]] == 1 &, Union[vars[[All, 1]]]] Positivity constraints:

varConstraints = Map[0 <= # <= 1 &, vars] Find a solution:

sol = Maximize[Join[{Total[vars]}, eqs, uniqueConstraints, varConstraints], vars, Integers] Replace the corresponding permutations and values:

varsToOne = Select[sol[], #[] == 1 &][[All, 1]];
vsol = Map[(#[] /. indexToPermRules) -> (#[] /. indexToMsValRules) &,
varsToOne];
ColumnForm[vsol]

TableForm[Transpose@{permsToMultiSetVals[[All, 1]],
permsToMultiSetVals[[All, 2]], permsToMultiSetVals[[All, 1]] /. vsol},
TableHeadings -> {None, {"permutations", "multiset values",
"solution replacements"}}, TableDepth -> 2] • A 0-1 ILP formulation does seem quite reasonable; that gets my upvote. – Daniel Lichtblau Nov 18 '15 at 19:06
• This is really a nice approach! – Per Alexandersson Nov 18 '15 at 19:09
• @DanielLichtblau Thanks -- I largely decided to look into this problem because of your comments and proposed solution. – Anton Antonov Nov 18 '15 at 19:28
• @PerAlexandersson Thanks! I would appreciate any comments or directions to how the complete solution would look like. – Anton Antonov Nov 18 '15 at 19:28
• @PerAlexandersson It seemed to me that a partially solved Sudoku puzzle is a good candidate for a practical example for the problem. My solution was derived from that point of view. I need to think about finding an example for the more general case. – Anton Antonov Nov 19 '15 at 17:51

This is probably not maximal efficiency but I think at least the narrowing it does is correct. Maybe.

I represent each partial maps as a set of two lists, where the first gives the domain values and the second gives the set of possible range values, regarded as unordered e.g. not corresponding in order to the associated domain values.

The idea is to iterate over pairs of such partial map. Any nontrivial intersection of domains and ranges gives a "reduction" or narrowing. Also such an intersection implies there are two complementary pairs and these too are in a sense narrowed as compared to their prior state. The task is to keep doing this narrowing until the map pair set stabilizes.

narrowMaps[lists : {{_List, _List} ..}] := Module[
{done = False, curlists = lists, newlists, l1, l2, vals1, vals2,
il1l2},
While[! done,
done = True;
newlists = Reap[
Do[
{l1, vals1} = curlists[[j]];
Do[
{l2, vals2} = curlists[[k]];
If[(il1l2 = Intersection[l1, l2]) === {}, Continue[]];
done = False;
Sow[{il1l2, Intersection[vals1, vals2]}];
curlists[[k]] = {Complement[l2, il1l2],
Complement[vals2, vals1]};
{l1, vals1} = {Complement[l1, il1l2],
Complement[vals1, vals2]};
, {k, j + 1, Length[curlists]}];
Sow[{l1, vals1}];
, {j, Length[curlists]}];
][[2, 1]];
If[! done, curlists = newlists];
];
newlists
]


I'll edit if and when some example or examples are proposed. Right now I simply illustrate with the one example that was mentioned in the post. Recall our domain->range set pairs consisted of {123, 213, 231, 312} associated to {1, 2, 0, 2} and {123, 321} associated with {3,0}. In my notation this set is represented as below.

mappairs = {{{123,213,231,312}, {1,2,0,2}}, {{123,321}, {3,0}}};


The example is now readily computed.

narrowMaps[ms1]

(* Out= {{{123}, {0}}, {{213, 231, 312}, {1, 2}}, {{321}, {3}}} *)


So the narrowing has determined that 123 must map to 0, ergo 321 must map to 3, and {213,231,312} maps to {1,2}.

--- edit ---

I realized just after posting that, due to the way we support multiple use of range values, this is not correct as written. In particular we know intersections force certain subsets of range values, but what I coded for complements is not right. Before attempting to correct it I'll first ask: do we always start with complete range sets counting by multiplicity? That is, if a range contains two 0's can I assume that each corresponds to exactly one (distinct) domain element? If so, the all that needs be done is to rework handling for proper accounting of multiplicity (well, "all" might be used in an optimistic sense there, but I don't think it should be too much of an alteration). If more is needed then best to have full details before undertaking a fix.

--- end edit ---

• Hi Daniel! Yes, the typical input should be that each object gives exactly one value under the (unknown) map. – Per Alexandersson Nov 17 '15 at 2:26