Opened 10 years ago

Closed 10 years ago

#11607 closed enhancement (fixed)

read constraints from linear program

Reported by: john_perry Owned by: ncohen
Priority: major Milestone: sage-4.8
Component: linear programming Keywords: sd32
Cc: Merged in: sage-4.8.alpha1
Authors: John Perry Reviewers: Nathann Cohen
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by ncohen)

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:

Attachments (2)

trac_11607_read_constraints_from_lp.patch (3.1 KB) - added by john_perry 10 years ago.
trac_11607_reviewer.patch (6.1 KB) - added by ncohen 10 years ago.

Download all attachments as: .zip

Change History (26)

comment:1 Changed 10 years ago by john_perry

  • Authors set to john_perry
  • Description modified (diff)
  • Status changed from new to needs_review

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

comment:2 Changed 10 years ago by john_perry

  • Keywords sd32 added

comment:3 follow-up: Changed 10 years ago by ncohen

  • Status changed from needs_review to needs_info

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

comment:4 in reply to: ↑ 3 ; follow-up: Changed 10 years ago by john_perry

Replying to ncohen:

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.

comment:5 in reply to: ↑ 4 ; follow-up: Changed 10 years ago by ncohen

Replying to john_perry:

Replying to ncohen: 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.

comment:6 in reply to: ↑ 5 ; follow-up: Changed 10 years ago by john_perry

Replying to ncohen:

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... :-)

comment:7 in reply to: ↑ 6 Changed 10 years ago by ncohen

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

comment:8 Changed 10 years ago by john_perry

  • Status changed from needs_info to needs_review

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().

comment:9 Changed 10 years ago by ncohen

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

comment:10 Changed 10 years ago by ncohen

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 ? :-)

comment:11 Changed 10 years ago by ncohen

  • Status changed from needs_review to needs_info

Changed 10 years ago by john_perry

comment:12 Changed 10 years ago by john_perry

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. :-)

comment:13 Changed 10 years ago by ncohen

  • Authors changed from john_perry to John Perry
  • Reviewers set to Nathann Cohen
  • Status changed from needs_info to needs_review

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:

Nathann

comment:14 Changed 10 years ago by ncohen

  • Description modified (diff)

comment:15 Changed 10 years ago by john_perry

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

comment:16 follow-up: Changed 10 years ago by ncohen

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.

comment:17 Changed 10 years ago by john_perry

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.

comment:18 in reply to: ↑ 16 Changed 10 years ago by john_perry

Replying to ncohen:

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.

comment:19 Changed 10 years ago by ncohen

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

Nathann

Changed 10 years ago by ncohen

comment:20 Changed 10 years ago by john_perry

  • Status changed from needs_review to positive_review

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...

comment:21 Changed 10 years ago by ncohen

Nice ! :-)

comment:22 Changed 10 years ago by jdemeyer

  • Milestone changed from sage-4.7.2 to sage-4.7.3

comment:23 Changed 10 years ago by jdemeyer

  • Milestone sage-4.7.3 deleted

Milestone sage-4.7.3 deleted

comment:24 Changed 10 years ago by jdemeyer

  • Merged in set to sage-4.8.alpha1
  • Milestone set to sage-4.8
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.