read constraints from linear program

[johnperry-math]

Sometimes a user might want to read & inspect the constraints in a linear program. For example, some code of mine generates a lot of linear programs (I described the project to Nathann once), and it would be nice if I could read the constraints off one program so as to combine it with another. Example usage would be:

sage: lp = MixedIntegerLinearProgram()
sage: lp.add_constraint(lp[0]-lp[1],min=1)
sage: lp.add_constraint(2*lp[1]-lp[0]-lp[2],max=1)
sage: lp.get_constraints()
[(1.0, ([1, 0],[-1.0, 1.0]), None),(None, ([2, 1, 0], [-1.0, 2.0, -1.0]), 1.0)]

Apply:

• attachment: trac_11607_read_constraints_from_lp.patch
• attachment: trac_11607_reviewer.patch

Component: linear programming

Keywords: sd32

Author: John Perry

Reviewer: Nathann Cohen

Merged: sage-4.8.alpha1

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

[johnperry-math]

comment:1

I've attached a patch that gives something resembling the desired functionality, as well as a couple of related, useful methods. I also changed the description, to mirror the desired result. (I do not really want to create a new symbolic expression describing a constraint at the Cython level.)

[johnperry-math]

Author: john_perry

[johnperry-math]

Changed keywords from none to sd32

[nathanncohen]

comment:3

Hello John !

I don't think the row and row_bounds method as they are implemented can be very useful to the user. They gave you access to the information of a constraint and take as an input the line number, but the user has for the moment no way whatsoever to know the number of a constraint. I mean, when would the user want to obtain the constraint corresponding to a constraint he does not know already, for instance obtained through your get_constraints method ?
I think those method, if you think they may beuseful by themselves, would make more sense in the generic_backend module. As all the other backends extend it, it would be available for ach backend, and could be used by your get_constraints method, it would make more sese at this level.

I have a similar remark for the get_constraints method. It returns a constraint as a list of indices and coefficient, but the indices have no meaning for the user who only deals with symbolic variables. Instead of returning vectors of indces and coefficient, which is the canonical form used by the solver's C API, what about returning a dictionary associating to each symbolic its coefficient ? This way the user can easily test whether the coefficient of one given variable is nonzero in the constraint, etc... :-)

Nathann

[johnperry-math]

comment:4

Replying to @nathanncohen:

I mean, when would the user want to obtain the constraint corresponding to a constraint he does not know already, for instance obtained through your get_constraints method ?

In my case, I have a program that generates constraints based on a polynomial system. I don't personally know what the constraints are, and in fact I want to read the generated constraints. For example, I might want to check in some easy circumstances that the program is generating the correct constraints.

I think those method, if you think they may beuseful by themselves, would make more sense in the generic_backend module.

Would that be immediately accessible to someone who has a MixedIntegerLinearProgram object?

Instead of returning vectors of indces and coefficient, which is the canonical form used by the solver's C API, what about returning a dictionary associating to each symbolic its coefficient ?

I didn't necessarily want the overhead of a dictionary. I could see adding that functionality, but if I were running a program where efficiency might be of some concern, generating a dictionary every time I do this seems a bit much.

Does that seem reasonable to you, or do you think I should generate the dictionary anyway? Efficiency here isn't an actual concern of mine at the moment, and might never be.

[nathanncohen]

comment:5

Replying to @johnperry-math:

Replying to @nathanncohen:
In my case, I have a program that generates constraints based on a polynomial system. I don't personally know what the constraints are, and in fact I want to read the generated constraints. For example, I might want to check in some easy circumstances that the program is generating the correct constraints.

I agree that in this case the MILP class should give you an easy way to check that, and to this end your get_constraints() method is a perfectly good answer. I still feel that the row(int index) and row_bounds(int index) feel more like backend methods, as they require an information that the user does not have : the indices of the rows. Even the number of constraints is not available right now :-)

Would that be immediately accessible to someone who has a MixedIntegerLinearProgram object?

No, in this case it would not, but this is my point exactly : should such functions be exposed, and if so in which aim ?

I didn't necessarily want the overhead of a dictionary. I could see adding that functionality, but if I were running a program where efficiency might be of some concern, generating a dictionary every time I do this seems a bit much.

When you talk of overhead, do you have memory or time in mind ? I do not think the memory cost would be that important, but as you repeatedly want to check for the existence of a given constraint, the "Set" structure which yields a log(n) check for existence seems more fitting than linearly exploring the cases each time a constraint is added.

Is it the idea of having a new element (the dictionary) attached to the MILP structure that you find ugly ? I have to admit it's not that sexy... To be perfectly frank I do not think it should be the MILP class' job to check whether the given set of constraints could be simplified, my way of doing it would be to build myself the set of constraints, which I would then use to build the LP (this would require the definition of __hash__ functions for several MILP-related classes). But then I try to keep in mind that we have different applications in mind, and that it can make much better sense in your situation :-)

If you think this would be a nice feature of the MILP class, however, I do think the dictionary trick is the best way to write it. :-)

Nathann

Does that seem reasonable to you, or do you think I should generate the dictionary anyway? Efficiency here isn't an actual concern of mine at the moment, and might never be.

[johnperry-math]

comment:6

Replying to @nathanncohen:

I agree that in this case the MILP class should give you an easy way to check that, and to this end your get_constraints() method is a perfectly good answer.

Yay :-)

I still feel that the row(int index) and row_bounds(int index) feel more like backend methods, as they require an information that the user does not have : the indices of the rows. Even the number of constraints is not available right now :-)

Actually, such a program might very well have counted the number of constraints added to the program. :-)

In any case, we can add a method to obtain this information, say, number_of_constraints(). We could go beyond this, and add to get_constraints() an optional argument, which_constraints, which is a list of constraints the user desires, rather than the entire list.

This doesn't feel backend-ish to me.

Would that be immediately accessible to someone who has a MixedIntegerLinearProgram object?

No, in this case it would not, but this is my point exactly : should such functions be exposed, and if so in which aim ?

I don't see why they shouldn't be exposed. As for the aim, IMHO it's always better to expose existing functionality that the user can use, at least if the functionality is read-only.

In my particular case, I might need to extend the constraints of one program by the constraints of another. This is probably not the only way to implement what I need done, but it is one way. Currently, I do it by building for each potential program a list of libsingular polynomials that correspond to each constraint, discard the programs after testing feasibility, then copying the equations of the best program into another one, using the coeff command. It's a fabulous bottleneck. The problem I'm working on generates a lot of programs.

(Notice I'm contradicting myself on efficiency -- I had forgotten how I might use these methods.)

I didn't necessarily want the overhead of a dictionary. I could see adding that functionality, but if I were running a program where efficiency might be of some concern, generating a dictionary every time I do this seems a bit much.

When you talk of overhead, do you have memory or time in mind ?

Either, but mostly time.

Is it the idea of having a new element (the dictionary) attached to the MILP structure that you find ugly ?

Here I think you're thinking I want this for #11606; that's not the case.

...my way of doing it would be to build myself the set of constraints, which I would then use to build the LP (this would require the definition of __hash__ functions for several MILP-related classes). But then I try to keep in mind that we have different applications in mind, and that it can make much better sense in your situation :-)

I'm perfectly open to any suggestions you might have; maybe we should take the discussion of the actual problem I'm working on into private email? I can send you a copy of the relevant sage code, but I don't know if you want to look at this monstrosity... :-)

[nathanncohen]

comment:7

Helloooooo !!!

In any case, we can add a method to obtain this information, say, number_of_constraints(). We could go beyond this, and add to get_constraints() an optional argument, which_constraints, which is a list of constraints the user desires, rather than the entire list.

I think "The Sage Way" would rather be something like "constraints()" instead of "get_constraints()" and an optional parameter, n for example, so that constraints(n = 10) returns constraint 10.

This doesn't feel backend-ish to me.

Yeah, well. Sometimes, takes me one week to notice I've been stupid. You're perfectly right.

I don't see why they shouldn't be exposed. As for the aim, IMHO it's always better to expose existing functionality that the user can use, at least if the functionality is read-only.

Right again. In my #687&%$&#687 of a head, they were exposed and useful in the GenericBackend class. Which was stupid, there's no way it would hurt to have them around in Python too.

In my particular case, I might need to extend the constraints of one program by the constraints of another. This is probably not the only way to implement what I need done, but it is one way. Currently, I do it by building for each potential program a list of libsingular polynomials that correspond to each constraint, discard the programs after testing feasibility, then copying the equations of the best program into another one, using the coeff command. It's a fabulous bottleneck. The problem I'm working on generates a lot of programs.

(Notice I'm contradicting myself on efficiency -- I had forgotten how I might use these methods.)

I hope we'll find a way to make this MILP class more useful to deal with such problems. But you will have to explain them to me, your use cases are very different :-)

Here I think you're thinking I want this for #11606; that's not the case.

I'm perfectly open to any suggestions you might have; maybe we should take the discussion of the actual problem I'm working on into private email? I can send you a copy of the relevant sage code, but I don't know if you want to look at this monstrosity... :-)

Why not ? I will not need to understand it all, but just the part in which the LP is generated. And "understanding" the LP isn't needed either, just how it is built :-)

Nathann

[johnperry-math]

comment:8

This version of the file changes the name of get_constraints() to constraints(). It adds an optional argument, indices, which defaults to None; if its value is None, then it returns a list of all constraints; otherwise, it returns a list of constraints indexed by None.

It also adds another function, number_of_constraints().

[nathanncohen]

comment:9

Hello again John !!!

I'm sorry for the time it takes to merge these patches, I'll be defending my PhD next month and I am crushed under the weight of the "urgent things" that should have been done already.

I read your patch again, and wondered at the terminology : "contraints" is used at some places, "row" at others. because I often mistake a row for a column I double-checked, but anyway your constraints method uses the row methods, so they are indeed returning the information on the constraints only, and not the variables as I had thought for a minute :-)

So, the thing is : what about having just one method constraints to do all the work ?

• constraints() would return all the constraints as a list, as it does now
• constraints(indices = [1,2,3]) would return constraints 1,2,3 as a list
• constraints(indices = 1) would return only one constraint, as the triples you used before. This way, all the job is done by one function, we stick to the terminology "constraints", and all is well under the sun ? :-)

It is one of the nice features of Python : having no types let us play a bit :-)

Tell me what you think of it. I'll try to make my next reviews muuuuch shorter !

Nathann

[nathanncohen]

comment:10

Ooops :

• constraints() would return all the constraints as a list, as it does now

• constraints(indices = [1,2,3]) would return constraints 1,2,3 as a list

• constraints(indices = 1) would return only one constraint, as the triples you used before. This way, all the job is done by one function, we stick to the terminology "constraints", and all is well under the sun ? :-)

[johnperry-math]

Attachment: trac_11607_read_constraints_from_lp.patch.gz

[johnperry-math]

comment:12

Sorry for the delay. I was working on some other things, including (but not limited to) that cvxopt change you mentioned in email.

Okay. I have made the changes you suggested. I went ahead and deleted the row() and row_bounds() methods, opting to stick with constraints(). On the other hand, I preserved the number_of_constraints() method, since that seems more straightforward.

One reason I had the row() and row_bounds() methods was for efficiency. But, you make a good point about row and constraint, and changing their names to constraint() and constraint_bounds() seemed likely to confuse the reader, for what is doubtless a microscopic gain in efficiency..

I'll have to revise #11606 based on this new patch, but don't forget it, either -- you haven't made any comments there in a while. :-)

[nathanncohen]

comment:13

Helloooooooooooo !!!

Sorry for all this time spent on that ticket. Your patch is nice, though it required a bit of other modifications:

• The documentation mentionned a former n parameter
• I reformatted this a bit to fit with the other methods
• I thought it was more natural to return only one triple and not a list of only one triple when indices is an integer
• There was something bad with the doctests : the order in which the values are returned (the indices and the coefficient) depends on the solver, so running tests on your file with CBC or CPLEX installed returned new errors. For this reason, I rewrote a bit the documentation so that the examples are at first not tested. In a second section, I wrote some code reordering the output before the tests so that the functions are indeed tested, and so that the tests work for any solver.
• While reviewing the code I found something really scary in CPLEX's interface : the constraints of the shape add_constraint(lb <= function <= ub) were actually totally wrong ! I added the line which was missing. This bug ony happened when adding both bounds at the same time (which you do in your example), so I hope no one was hurt :-/. Anyway it is far better with this fixed !

I attach to this ticket my reviewer's patch. If you agree with the changes, you can set this ticket to "positive review" :-)

Apply:

• attachment: trac_11607_read_constraints_from_lp.patch
• attachment: trac_11607_reviewer.patch

Nathann

[nathanncohen]

Reviewer: Nathann Cohen

[nathanncohen]

Changed author from john_perry to John Perry

[johnperry-math]

comment:15

Should I be concerned about the lines in the doctests of your patch that say, "not tested"? For example,

    To obtain the list of all constraints:: 

            sage: p.constraints()          # not tested

That's my only concern with your patch (I haven't tested it yet though -- will do so in a second).

[nathanncohen]

comment:16

I explained why in my description of the patch. These lines are indeed not tested, because that would create wrong errors when using different solvers, but the returned values are tested for correctness in a TESTS section later in the patch.

[johnperry-math]

comment:17

A problem: line 457,

sage: p.constraints(0)

should actually read,

sage: reorder_constraint(p.constraints(0))

If I make that one change, doctests pass.

[johnperry-math]

comment:18

Replying to @nathanncohen:

I explained why in my description of the patch. These lines are indeed not tested, because that would create wrong errors when using different solvers, but the returned values are tested for correctness in a TESTS section later in the patch.

You're right. Sorry about that.

[nathanncohen]

comment:19

Oopsssssss !! You're right too ! patch updated :-)

Nathann

[johnperry-math]

comment:20

Attachment: trac_11607_reviewer.patch.gz

I'm giving it a positive review, based on what Nathann wrote above:

If you agree with the changes, you can set this ticket to "positive review" :-)

Doctests pass...

[nathanncohen]

comment:21

Nice ! :-)

[jdemeyer]

Milestone sage-4.7.3 deleted

[jdemeyer]

Merged: sage-4.8.alpha1