Originally posted by: Jernej1

Hello!

I am working on a binary ILP with some very simple constraints. I am generating the LP files with SageMath (GLPK?) and solving them with CPLEX. 

I have recently started porting the problem generation part directly into CPLEX's C++ interface and exporting the LP files with IloCplex::exportModel. 

To my surprise, it takes significantly longer to solve problems that were exported with CPLEX, to the point of suspecting there is some other issue with the way I am generating the problems.

I've created two minimal examples (L1min.lp, L2min.lp created with CPLEX and SageMath respectively) where one can already see the difference in the solution time. The run time is indeed insignificant on them, but it shows the general structure of the problem I am dealing with.

The two larger instances (L1.lp, L2.lp) already offer a significant difference in execution time.

I am aware that this phenomenon is to be expected (different variable orderings giving different execution times) but the problem is that I was not able to find an example where generating models with CPLEX would yield an improved run time.

Before I resign into the performance setback I would at least like to make sure I am perhaps not mismodeling the problems.

Can someone offer some insight into this issue? Is there a way to quickly verify that L1.lp,L2.lp are indeed isomorphic and offer some guideline how to achieve a similar performance?

Given that these are binary ILP's with binary coefficients I could probably use graph theory to figure out isomorphism, but I am happy enough if someone could give a quick glance and sense if something smells fishy.

How sure are you that these models are the same?

The first two constraints in L1.lp are

 c1:    x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 >= 1
 c2:    x1 + x2 + x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 >= 1

The first two constraints in L2.lp are

 r_1: - x_11 - x_10 - x_9 - x_8 - x_7 - x_6 - x_5 - x_4 - x_3 - x_2 - x_1 <= -1
 r_2: - x_20 - x_19 - x_18 - x_17 - x_16 - x_15 - x_14 - x_13 - x_12 - x_6 - x_1 <= -1

r_2 references x_6 while c2 references x2, the rest of the rows are the same. Of course, this may  just be due to permutation of variables and the respective blindly naming variables x1,...,xn. However, before digging deeper, I'd rather double-check. Do you have a way to create unique names for the variables so that you can be sure variables with same names in the different models mean the same thing?

Originally posted by: Jernej1

Daniel,

thank you for taking the time to check the LP's. Given your suggestion I'll try to make a 100% definite isomorphism check (good point with trying to name the variables equivalently) and come back once this is fully settled.

Originally posted by: Jernej1

Hello,

I just assured myself that the two problems are isomorphic and is just an issue of how variables are being named. The variables in the LP generated by SageMath are indexed by whatever the underlying indexation of Pythonic dicts is, and hence I don't really expect to be able to reproduce that in the coresponding C++ program.

The runtime difference between L1.lp and L2.lp is of magnitude 10x ticks. 

What ways could I explore to solve L1.lp faster?

I am observing strange things here:

I randomly permuted rows and columns of both models 30 times and solved the resulting models. For one model I get rather fast consistent running times. For the other model the running times very a lot, between 30 and thousands of seconds. That still makes me wonder whether the two models are indeed the same.

Can you try the following: Currently your model does not have any objective function. Can you add as objective

  minimize sum (j+1) * xj

to both models and export/solve them again. That should also make it very simple to check for isomorphism since the column permutation is obvious and you only have to check whether there is a row-permutation that transform one model into the other. Also, with this artificial objective L1 solved rather quickly here.

Originally posted by: Jernej1

Thank you for taking the effort. I also observed a similar sampling inconsistency with solution times and hence raised this issue here. By forcing a different ordering of the constraints being added to the Sagemath model I convinced myself that I get exactly the same LP files. I am not any mroe.

Adding the objective function makes the solution time of L2.lp  much shorter, but L1.lp still seems to be struggling a lot.

Unfortunately I am not sure I see the obvious isomorphism that the objective function brings in. Could you help me a bit with this? Do you happen to see it, how can I produce the bijection?

Thank you,

Jernej

The idea is this: Each variable is uniquely identified by its objective coefficient. Consider for example the variables with objective coefficient 1. Each model has exactly one such variable. Let j1a denote that variable in the first and j1b denote that variable in the second model. For the two models to be isomorphic, we have to map j1a to j1b. There is no other choice. So we can read the required mapping directly from the objective coefficients. Now, in both models we can rename the variables so that the variable with objective coefficient k is called "jk".

Now the variables in the two models have the same names and variables with same names are mapped to each other. To prove or disprove isomorphism of the two models, we just have to find for each constraint in the first model the same constraint in the second model.

I did all the above using the attached Python script. And (modulo bugs in my script), the two models are not somorphic. The '>=' constraints in the two models are the same, but the '=' constraints seem to be different. For example, both models have exactly one equality constraint with right-hand side 8. In the first model, this constraint reads:

  c1029: j5 j13 j21 j29 j37 j45 j53 j61 j69 j77 j85 j93 j101 j109 j117 j125 j133 j141 j149 j157 j165 j173 j181 j189 j197 j205 j213 j221 j229 j237 j245 j253 j261 j269 j277 j285 j293 j301 j309 j317 j325 j333 j341 j349 j357 j365 j373 j381 j389 j397 j405 j413 j421 j429 j437 j445 j453 j461 j469 j477 j485 j493 j501 j509 j517 j525 j533 j541 j549 j557 j565 j573 j581 j589 j597 j605 j613 j621 j629 j637 j645 j653 j661 j669 j677 j685 j693 j701 j709 j717 j725 j733 j741 j749 j757 j765 j773 j781 j789 j797 j805 j813 j821 j829 j837 j845 j853 j861 j869 j877 j885 j893 j901 j909 j917 j925 j933 j941 j949 j957 j965 j973 j981 j989 j997 j1005 j1013 j1021 = 8

While in the second model the constraint reads:

  c1031: j513 j514 j515 j516 j517 j518 j519 j520 j521 j522 j523 j524 j525 j526 j527 j528 j529 j530 j531 j532 j533 j534 j535 j536 j537 j538 j539 j540 j541 j542 j543 j544 j545 j546 j547 j548 j549 j550 j551 j552 j553 j554 j555 j556 j557 j558 j559 j560 j561 j562 j563 j564 j565 j566 j567 j568 j569 j570 j571 j572 j573 j574 j575 j576 j577 j578 j579 j580 j581 j582 j583 j584 j585 j586 j587 j588 j589 j590 j591 j592 j593 j594 j595 j596 j597 j598 j599 j600 j601 j602 j603 j604 j605 j606 j607 j608 j609 j610 j611 j612 j613 j614 j615 j616 j617 j618 j619 j620 j621 j622 j623 j624 j625 j626 j627 j628 j629 j630 j631 j632 j633 j634 j635 j636 j637 j638 j639 j640 = 8

As you can see, the involved variables are completely different.

You can even spot the difference in the original models. In L2.lp the constraint is

x2 + ... = 8

and the objective coefficient of x2 is 513.

In L1.lp the variable with objective coefficient 513 is x11 and that variable does not appear in the =8 constraint at all.

With that I am pretty sure that the two models are not permutations of each other.

Originally posted by: Jernej1

Daniel,

thank you very much for such an extensive effort.

I am really baffled by this as I was sure the models are isomorphic. I tested both SageMath/C++ generators on smaller instances and made sure the same variables are being used while the model is built. And everything works out well for small models and the isomorphism is also very much apparent there (speaking of problems with 4 and 8 variables).

At this point I can only assume there is a weird bug somewhere in how the SageMath model is being written in the LP file or perhaps I am generating the problem wrong in C++. Since it is simple to understand I will post it here and hopefully someone can confirm that it works as desired.

The problem is a variant of the binary covering code problem. You are given the set of all binary strings of length 10 and you want to find a subset S of strings so that every other binary string differs in at most one coordinate from some element in S. This is what the first 1024 constraints are for.

For the second set of constraints, you are given a list of x_0,x_1,...,x_{2^K} integers for example 42,9,9,13,8,13,13,12.

The condition dictates that for any binary K-string F the number of elements of S that have F as a prefix is precisely x_F, where I am treating F as both a string and the respective integer representation here. These are the remaining 8 conditions in the attached LP file. The C++ program that generates this problem in CPLEX is

#define N 10
#define K 3 
#define CONF_SIZE ( 1<<K )

IloEnv env;

/* This function returns 0 if it could not prove that the given
 * configuration is infeasible, nonzero otherwise
 */

static unsigned isInfeasible(const unsigned d[CONF_SIZE]) {

    unsigned i,j;        
    int status;

    try {

        IloModel model(env);
        IloNumVarArray var(env);
        IloRangeArray con(env);
        IloCplex cplex(model);

        for (i = 0; i < (1U<<N) ; i++) {
            var.add( IloBoolVar(env)); 
        }

        /* Covering conditions */        
        for (i = 0; i < (1U << N) ; i++) {
            IloExpr expr(env);
            expr =  var[i];
            for (j = 0; j < N; j++) {
                expr += var[i ^ (1<<j)];
            }
            con.add( expr >= 1 );         
        }

        /* Conditions related to the specific configuration */
        for (i = 0; i < CONF_SIZE; i++) {
            IloExpr expr(env);
            for (j = 0; j < (1U << (N-K) ); j++) {
                expr += var[ (j << K) | i ];
            }
            con.add( expr == d[i] );
        }

        model.add(con);
        cplex.solve();

        if ( cplex.getStatus() == IloAlgorithm::Infeasible ) {
            return 1;
        }
    } catch(IloException& e) {
        std::cerr << "Concert exception caught: " << e << std::endl;
    }

    catch(...) {
        std::cerr << "Unknown exception caught" << std::endl;
    }
    return 0;
}

Can you confirm that the following indeed matches the description of the model?

Thank you again for your assistance!

As far as I can tell, your C++ is correct.

I don't know SageMath/C++ but if the code is as short as your C++ code, could you post that as well? Maybe we can spot the difference.

Maybe also post the C++ and SageMath code that creates this fake objective function.

Originally posted by: Jernej1

Hey there,

thanks for confirming. Yes I did not want to post SageMath code as I thought it might be too off topic for this forum. Following is the code. Few remarks

- Whether we use the equality constraint on prefixes/suffixes should be symmetric

- I've checked that the generated constraint use the same numbers (binary strings) as the C++ program

- The function accepts a dictionary mapping K-binary strings to the respective integer values as described above.

N = 10

def writeLP(d):
    # SageMath CubeGraph gives us the following two properties 
    # for v in G.vertices() - will iterate over all N-binary strings 
    # for a N-binary string f, G[f] gives a list of its neighbors

    G = graphs.CubeGraph(N)
    p=MixedIntegerLinearProgram(maximization=False, solver='CPLEX')
    b=p.new_variable(binary=True)

    for v in G.vertices():
        p.add_constraint(b[v]+p.sum([b[u] for u in G[v]]) >= 1)

    for key in a:
        p.add_constraint ( sum(b[v] for v in G if v[:len(key)] == key) == d[key]) 

    p.write_lp('-'.join(str(d[key]) for key in sorted(d.keys()))+ '.lp')

I wonder whether there may be a problem with the translation of the right-hand sides of the equality constraints. As far as I can tell, in your SageMath code you read binary strings from left to right, while in the C++ code you definitely read the strings from right to left.

That means that what is d["100"] in SageMath, must be at d[001] in the C++ code, correct? Can you double check that?

Also, the idea of marking the variables by unique objective coefficients may need to be tweaked since the binary strings are read in different directions in the two models. Hence in the C++ model, the coefficient for a variable should be derived from the value of the reversed binary string.

Originally posted by: Jernej1

Daniel,

thanks for all the effort!

As I mentioned in the above post the code indeed use different "endianess". But this in itself is not a problem. If S is a covering code, then so is the code obtained by reversing the strings in it. So if you set your condition for prefixes or suffixes is irrelevant.

What is relevant (and I completely missed that one) is that if you set the condition on the prefixes (instead of suffixes) you also need to reverse the indexes for the elements in d. For the above example, the element at position 1 is actually at position 8, an element at position 6 is at position 3 and so on...

So I believe this to be the actual issue.

Thanks again for all the effort in resolving this. 

The following C++ program produces a model that is isomorphic to the latest L2.lp that you posted. Maybe you can take it from there to figure out what your C++ code does differently and whether you have to update your C++ code or the SageMath code. I  would put my money on the order of the elements in the 'd' array (which is called 'lim' in my code.

#include <ilcplex/ilocplex.h>

#define N 10
#define K 3 
#define CONF_SIZE ( 1<<K )

using std::cout;
using std::endl;

unsigned reversebits(unsigned i) {
   unsigned res = 0U;
   for (unsigned j = 0, mask = 1U << (N - 1); j < N; ++j, mask /= 2)
      res |= ((i & (1U << j)) >> j) * mask;
   return res;
}

int
main(void)
{
   static int const lim[] = { 42, 8, 9, 13, 9, 13, 13, 12 };

   IloEnv env;
   IloModel model(env);
   
   IloNumVarArray x(env);
   IloExpr obj(env);
   for (unsigned i = 0; i < (1U << N); ++i) {
      IloBoolVar v(env);
      x.add(v);
      obj += int(1 + reversebits(i)) * v;
   }
   model.add(IloMinimize(env, obj));

   for (unsigned i = 0; i < (1U << N); ++i) {
      IloExpr sum = x[i];
      for (unsigned j = 0; j < N; ++j)
         sum += x[i ^ (1U << j)];
      model.add(sum >= 1);
   }

   for (unsigned i = 0; i < CONF_SIZE; ++i) { // All binary strings of length K
      IloExpr sum(env);
      for (unsigned j = 0; j < (1U << (N-K) ); ++j) {
         sum += x[(j << K) | i];
      }
      model.add(sum == lim[i]);
   }

   IloCplex cplex(model);
   cplex.exportModel("model.lp");
   env.end();

   return 0;
}