a bijectionist's toolkit

[mantepse]

We provide a toolkit for the combinatorialist to help find functions ("statistics") s: A -> Z and bijections A -> B given sequences of finite sets A and B that satisfy various constraints.

Depends on #34881 Depends on #34878

CC: @stumpc5 @tscrim @fchapoton

Component: combinatorics

Author: Alexander Grosz, Tobias Kietreiber, Stephan Pfannerer, Martin Rubey

Branch/Commit: u/mantepse/a_bijectionist_s_toolkit @ 53506aa

Reviewer: Matthias Koeppe, ...

Issue created by migration from https://trac.sagemath.org/ticket/33238

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

db850f0

add and remove constraints instead of copying the whole program

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

ec3271c

fix problem with SCIP, add timings

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from db850f0 to ec3271c

[mkoeppe]

comment:37

+            d = minimal_subdistribution.get_values(D)  # a dict from P to multiplicities
[...]
+                support = [X[p] for p in P if d[p]]

I would recommend to use the recently added convert parameter of get_values https://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html#sage.numerical.mip.MixedIntegerLinearProgram.get_values, which makes such tests much more robust.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from ec3271c to dfdc6ba

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

dfdc6ba

use convert for MixedIntegerLinearProgram.get_values

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

d0836cb

include information on MILP also in public documentation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from dfdc6ba to d0836cb

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

563af98

remove unnecessary copy

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from d0836cb to 563af98

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

590bae2

add possibility to constrain to involutions

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 563af98 to 590bae2

[mantepse]

comment:42

See https://mathoverflow.net/questions/437261/standard-or-good-names-for-relations-between-maps for a question on terminology.

If you have a nice test for the constraint s o psi o s = phi, with non-identity psi and phi, I'd love to add it.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

da0655b

correct indentation

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 590bae2 to da0655b

[mantepse]

comment:44

A question about the general design: currently there is a helper class _BijectionistMILP. I am not quite sure whether it might be better to merge it into the main class or not.

The idea in creating it was to separate the MILP from the user interface. More precisely, the methods Bijectionist.set_* in the main class are responsible to normalizing and organizing the input, whereas the methods _BijectionistMILP.add_*_constraint in the helper class actually turn the input into constraints.

Also, the helper class organizes a cache of solutions, and, whenever a new solution is found, a veto constraint is added to the MILP. Put differently, the set of solutions in the cache together with the set of solutions in the MILP is always constant.

I think the current solution is OK, but I am open to improvements.

[mantepse]

comment:45

Concerning the foreign_latex plugin: it is triggered by overset because of the regexp over[^l]. I believe that my use of overset is valid, I do not know of a better construct here. In any case, the commutative diagram is typpeset correctly.

[mantepse]

comment:46

Possibly, using the indices argument of new_variable, I could slightly simplify parts of the code. Currently, variables are added as needed, and blocks where this is not allowed are guarded by checking whether the number of variables changed.

I am not sure how easy it is to determine the allowed indices - they are a subset of P x Z, where P is the set of constant blocks and Z is the set of values.

[mkoeppe]

comment:47

Replying to Martin Rubey:

A question about the general design: currently there is a helper class _BijectionistMILP. I am not quite sure whether it might be better to merge it into the main class or not.

It may help to think about the overall structure of the code by considering whether 1. other enumeration technologies could be supported, or 2. whether some of your code could be generalized.

1. Alternative techniques for the 0/1 enumeration: lattice point enumeration (Polyhedron.integral_points()), SAT (e.g. #33461 cvc5).

2. Here you enumerate 0/1 solutions by repeatedly solving a MIP with no-good cuts added. If you generalize the form of your cuts a bit (so that their validity does not depend on equal cardinality), then you'd have a general algorithm for enumeration that could take as input any 0/1 IP -- which could be a useful addition to the Sage library. (It could later be improved by solver-specific techniques such as accessing the solver's solution pools, or using solver callbacks.)

[mkoeppe]

comment:48

Replying to Martin Rubey:

Concerning the foreign_latex plugin: it is triggered by overset because of the regexp over[^l]. I believe that my use of overset is valid, I do not know of a better construct here. In any case, the commutative diagram is typpeset correctly.

Yes, it's safe to ignore this warning. When we switch to GitHub, the patchbot will be retired anyway. The checks on GH Actions (tox -e relint) have a better pattern for foreign_latex, see src/.relint.yml.

[mantepse]

comment:49

Replying to Matthias Köppe:

Replying to Martin Rubey:

A question about the general design: currently there is a helper class _BijectionistMILP. I am not quite sure whether it might be better to merge it into the main class or not.

It may help to think about the overall structure of the code by considering whether 1. other enumeration technologies could be supported, or 2. whether some of your code could be generalized.

1. Alternative techniques for the 0/1 enumeration: lattice point enumeration (Polyhedron.integral_points()), SAT (e.g. #33461 cvc5).

2. Here you enumerate 0/1 solutions by repeatedly solving a MIP with no-good cuts added. If you generalize the form of your cuts a bit (so that their validity does not depend on equal cardinality), then you'd have a general algorithm for enumeration that could take as input any 0/1 IP -- which could be a useful addition to the Sage library. (It could later be improved by solver-specific techniques such as accessing the solver's solution pools, or using solver callbacks.)

Could you please expand? I do not really understand what you mean. In particular, I do not know what a cut is, and what you mean with "equal cardinality". Is the idea that to have a class that gets you all solutions for a MIP, instead of just one?

I admit that I am slightly afraid of increasing the scope of the project at this stage. Did you have a look at the MILPs solved? (some are documented in solutions_iterator, there are some more in _forced_constant_blocks, minimal_subdistributions_iterator and minimal_subdistributions_block_iterator)

[mantepse]

comment:50

Please excuse, I have yet another question. I was trying to simplify evaluate_linear_function, is there code to check whether a given set of values satisfies a constraint?

Besides, I asked chatGPT. It answered as follows:

--- begin quote ---

Yes, you can use the MILP class from the MixedIntegerLinearProgram module in SageMath to check whether a given set of values satisfies a constraint of a MILP.

Here is an example of how you can use the MILP class to check whether a given set of values satisfies a constraint of a MILP:

from sage.numerical.mip import MixedIntegerLinearProgram

# Create a new MILP
milp = MixedIntegerLinearProgram()

# Set up variables
x = milp.new_variable(real=True, nonnegative=True)
y = milp.new_variable(real=True, nonnegative=True)

# Set up constraints
milp.add_constraint(x + y <= 1)
milp.add_constraint(x - y >= 0)

# Check whether the given set of values satisfies the constraints
satisfies_constraints = milp.check_feasibility([x: 0.5, y: 0.5])

print(satisfies_constraints)  # prints True

Could you show me the definition of check_feasibility?

Certainly! Here is the definition of the check_feasibility method of the MILP class from the MixedIntegerLinearProgram module in SageMath:

def check_feasibility(self, values=None, variable_order=None, solver=None):
    """
    Check the feasibility of a given set of values.

    INPUT:

    - ``values`` -- a dictionary of variable names to values.

    - ``variable_order`` -- a list of variable names in the order in which they
      should be presented to the solver. If not specified, the default
      order is used.

    - ``solver`` -- the name of the solver to use. If not specified, the
      default solver is used.

    OUTPUT:

    Returns True if the given set of values satisfies all of the constraints
    of the MILP, and False otherwise.

    EXAMPLES:

    Check the feasibility of a given set of values::

        sage: from sage.numerical.mip import MixedIntegerLinearProgram
        sage: milp = MixedIntegerLinearProgram()
        sage: x = milp.new_variable(real=True, nonnegative=True)
        sage: y = milp.new_variable(real=True, nonnegative=True)
        sage: milp.add_constraint(x + y <= 1)
        sage: milp.add_constraint(x - y >= 0)
        sage: milp.check_feasibility([x: 0.5, y: 0.5])
        True

    """
    # ...

--- end quote --- Unfortunately, this method does not exist in the current sagemath develop branch :-(

[mkoeppe]

comment:51

Replying to Martin Rubey:

Replying to Matthias Köppe:

2. Here you enumerate 0/1 solutions by repeatedly solving a MIP with no-good cuts added. If you generalize the form of your cuts a bit (so that their validity does not depend on equal cardinality), then you'd have a general algorithm for enumeration that could take as input any 0/1 IP

Could you please expand? I do not really understand what you mean. In particular, I do not know what a cut is, and what you mean with "equal cardinality".

Sorry for the jargon. "Cut" is short for "inequality". What you call a "veto inequality" is commonly called a "no-good cut".

When you have a solution x1=1, x2=1, x3=0, x4=0, you are generating an inequality x1+x2 <= 1 in order to cut off this solution, but this is only valid because you know that x1+x2+x3+x4 = 2. The more general form of the no-good cut would be (1-x1)+(1-x2)+x3+x4 >= 1, which is linearly equivalent to your inequality given x1+x2+x3+x4 = 2.

[mantepse]

comment:52

Is it really the case that LinearConstraint does not know its variables, other than by integer index, and that, given an index, there is no way to obtain the key?

Explicitly:

sage: m = MixedIntegerLinearProgram()
sage: x = m.new_variable(indices=["a", "b", "c"])
sage: c = x["a"] + 2*x["b"] >= 3
sage: low, f = next(c.inequalities())
sage: x.items()
dict_items([('a', x_0), ('b', x_1), ('c', x_2)])
sage: f.dict()
{0: 1.0, 1: 2.0}

What I'd like to do is to "evaluate" f, and I think that this should be a method in LinearConstraint, possibly it should be the subs method, but I am unsure, because I'd like it to return True or False if the constraint is satisfied by the given values.

[mkoeppe]

comment:53

Replying to Martin Rubey:

Is it really the case that LinearConstraint does not know its variables, other than by integer index, and that, given an index, there is no way to obtain the key?

Yes, unfortunately. That's basically #26302.

[mkoeppe]

comment:54

see also #20664.

[mantepse]

comment:55

I admit I am confused by the multitude of tickets, and I also lack understanding the basic design of the frontend. In particular, it is unclear to me, why some classes and methods seem to use integers to index the variables ("columns", I guess), such as LinearFunction.dict, and some use the indices by the user, such as get_values.

Wouldn't it be better to use the user supplied indices always, except in the backend?

[mkoeppe]

comment:56

Replying to Martin Rubey:

I also lack understanding the basic design of the frontend.

Well, it's inconsistent. In the long run, I hope that we can replace our frontend with something more disciplined. My current candidate for this is CVXPY (#34251).

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from da0655b to 03b1fb8

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

03b1fb8

better handling of empty constraints

[mkoeppe]

Dependencies: #34881

[mkoeppe]

comment:59

Replying to Martin Rubey:

A question about the general design: currently there is a helper class _BijectionistMILP. I am not quite sure whether it might be better to merge it into the main class or not.

Some more comments on this:

• the steps of _initialize_new_bmilp could be merged into _BijectionistMILP.__init__
• _show_bmilp would make more sense as a method of _BijectionistMILP (it breaks the abstraction of _BijectionistMILP)
• also _solution could be changed to a (public) method of _BijectionistMILP; and perhaps there could be a slightly more high-level function so that various Bijectionist methods do not have to repeat the "try ... solve ... except MIPSolverException ... else solution" pattern

[mkoeppe]

comment:60

I would also suggest to introduce a method in _BijectionistMILP that returns the variable _x[p, z] so that you can remove the direct use of this attribute from the Bijectionist class.

[mkoeppe]

Changed dependencies from #34881 to #34881, #34878

[mkoeppe]

comment:62

By the way, in the benchmark example (bij), instead of using list(bij.solutions_iterator()), the following works:

sage: bmilp = bij._initialize_new_bmilp()
sage: %time p = bmilp.milp.polyhedron(base_ring=QQ)
CPU times: user 126 ms, sys: 13.9 ms, total: 140 ms
Wall time: 164 ms
sage: p.is_lattice_polytope()
True
sage: len(p.vertices())
504

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

a3364de

move _show to the _BijectionistMILP class

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 03b1fb8 to a3364de

[mantepse]

comment:64

Replying to Matthias Köppe:

By the way, in the benchmark example (bij), instead of using list(bij.solutions_iterator()), the following works:

sage: bmilp = bij._initialize_new_bmilp()
sage: %time p = bmilp.milp.polyhedron(base_ring=QQ)
CPU times: user 126 ms, sys: 13.9 ms, total: 140 ms
Wall time: 164 ms
sage: p.is_lattice_polytope()
True
sage: len(p.vertices())
504

Wow, this is very impressive. Why is this so much faster? Can I always use this method instead?

[mantepse]

comment:65

I guess that, in general, we want integral_points, right?

Oh, but then we are back to a long time. Or do I misunderstand something? But that's strange, it should be fast to filter out non-integral points.

[mkoeppe]

comment:66

It's problem-dependent how well this works.

[mkoeppe]

comment:67

For a 0/1 problem, every 0/1 solution is a vertex solution of any polyhedral relaxation that keeps the constraints 0<=x_i<=1.

So you can enumerate the 0/1 solutions by filtering all vertex solutions to any such polyhedral relaxation. Whether this works well depends on the relaxation. In the example, your problem formulation happens to be perfect: Every vertex is a 0/1 solution. In general, you can have exponentially more vertex solutions than 0/1 solutions. In that case, this algorithm will not work so well.

[mkoeppe]

comment:68

Replying to Martin Rubey:

I guess that, in general, we want integral_points, right?

Oh, but then we are back to a long time.

The default implementation of integral_points is not very good for this use case. But after sage -i pynormaliz, you can do

sage: bmilp = bij._initialize_new_bmilp()
sage: %time p = bmilp.milp.polyhedron(base_ring=QQ, backend='normaliz')
CPU times: user 110 ms, sys: 3.31 ms, total: 113 ms
Wall time: 112 ms
sage: %time p = bmilp.milp.polyhedron(base_ring=QQ, backend='normaliz').integral_points()
CPU times: user 317 ms, sys: 35.3 ms, total: 352 ms
Wall time: 380 ms

In this example (perfect formulation), it only adds overhead.

It would be interesting to investigate bijectionist examples in which the IP solver actually has to do something nontrivial.

[mantepse]

comment:69

In any case I would like to structure the program so that we can try either method easily. Possibly this guides us towards a good design?

I should add, the most important use case is probably minimal_subdistributions_iterator. Will it work there, too?

Unfortunately, I am running out of time myself :-(

[mkoeppe]

comment:70

Replying to Martin Rubey:

In any case I would like to structure the program so that we can try either method easily. Possibly this guides us towards a good design?

Yes, I think my suggestions in comment:59, in particular avoiding use of MIPSolverException in the main class, will go in this direction

[mkoeppe]

comment:71

Replying to Martin Rubey:

the most important use case is probably minimal_subdistributions_iterator. Will it work there, too?

I haven't looked at the IP formulation that you use there, but I see that you have some general integer variables in addition to 0/1 variables. Is it true that a choice of the values of the 0/1 variables uniquely determines the values of the other variables? Then the same method applies.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

9d83a6c	remove unused code for backends that do not support removing constraints
864029d	make _solution a public method of _BijectionistMILP, simplify solve

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a3364de to 864029d

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 864029d to 64101a8

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

64101a8

eliminate _initialize_new_bmilp

[mantepse]

comment:74

I have now completed most of your suggestions. The remaining ones are:

• perhaps there could be a slightly more high-level function [replacing _BijectionistMILP.solution] so that various Bijectionist methods do not have to repeat the try ... solve ... except MIPSolverException ... else solution pattern

• I would also suggest to introduce a method in _BijectionistMILP that returns the variable _x[p, z] so that you can remove the direct use of this attribute from the Bijectionist class.

I am not completely sure how exactly I should go about these. I am still not convinced that it is a good idea to have the helper class at all. The distinction between "temporary constraints", which are built in the main class and passed to _BijectionistMILP.solve and "permanent constraints", which are built by _BijectionistMILP itself, doesn't look completely right to me.

Also, it seems a bit strange that _BijectionistMILP calls methods of the main class, _preprocess_intertwining_relations and _compute_possible_block_values. Maybe these should be moved to the helper class also?

On a quite different level, I am pretty sure that the logic of _forced_constant_blocks can be simplified. Apart from the main routine, note that at the very end we do self.set_constant_blocks, which invalidates the MILP and therefore all solutions in the cache.

Yet another thing, which is truly a detail: the method _is_solution really should be a method of LinearConstraint.

In summary, there is still a lot of work that could be done, but I am not sure whether the right time for this is now.