## linear programming

Librairie C++ de calcul formel/ C++ symbolic computation library

parisse
Messages : 5082
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: linear programming

Thank you, I'm incorporating your changes.
hevea complains about using constructs like _\mathrm{eq}, could you replace with_{\mathrm{eq}} for the next time (I removed \mathrm this time)?

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

I'll do that. Thanks for the quick upgrade!

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

I am preparing an (significantly) improved version of the code, which I will post in few days, when I test and debug it.
Bernard, what do you think about enabling lpsolve to use facilities from GLPK (GNU Linear Programming Kit) library when computing in floating point arithmetic? It would mean a larger dependency list, however, unless the compilation with or without GLPK is made optional. The alternative is to keep the current implementation of the primal-dual affine scaling method (which is improved, by the way, but still not far from a "toy" implementation). If the latter is preferable, would computing with c doubles instead of gens (containing doubles) make the algorithm faster?

parisse
Messages : 5082
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: linear programming

No problem to add a GPL library, especially if it's optionnal for better performances.
And yes, working with double instead of gen is significantly faster (and a little more precise).

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

Hello Bernard, I have a question. When specifying integer and binary variables in lpsolve when problem is given in matrix form, user has to enter indices of the corresponding variables. I want to make these indices 0-based for Xcas and Maple mode and 1-based for Mupad mode (as it is the default starting ordinal in these modes). How do I do that? I want this feature to be consistent with mode chosen.

parisse
Messages : 5082
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: linear programming

Did that change in Maple? In Maple V, indices were 1-based.
With the latest versions of giac, you can get the start index for arrays by calling array_start(contextptr)

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

My mistake, Maple indices are indeed 1-based. Thanks for pointing out the array_start function, it's just what I needed.

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

To use array_start, I have to install latest giac library. But apt-get in Ubuntu does not allow me to use your testing branch:
deb http://www-fourier.ujf-grenoble.fr/~parisse/debian/ testing main
It says that Release file is missing. Should I compile giac on my own instead? Or there's an easier way?

parisse
Messages : 5082
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: linear programming

Get the .deb and install it with dpkg -i
Or use the code of array_start

Code : Tout sélectionner

  int array_start(GIAC_CONTEXT){
if (contextptr && contextptr->globalptr){
return (!contextptr->globalptr->_python_compat_ && (contextptr->globalptr->_xcas_mode_ || absint(contextptr->globalptr->_calc_mode_)==38))?1:0;
}
return (!_python_compat_ && (_xcas_mode_ || absint(_calc_mode_)==38))?1:0;
}


lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

Here is the new version of the code. It was almost completely rewritten to incorporate better techniques and be more readable and easier to maintain. New features and improvements are listed below.

1. there are new options: lp_maxcuts, lp_gaptolerance, lp_nodeselect, lp_varselect, lp_iterationlimit, lp_timelimit and lp_verbose (see documentation entry for explanations). Mostly these are used to control branch&bound process. There are also new keywords used as values for some of these options. All new keywords (options and values) are listed at the top of lpsolve.h and they should be moved to dispatch.h. Some options became obsolete and should be removed. These are lp_integertolerance, lp_initialpoint and lp_variables. Bernard, could you please take care of this?
2. lpsolve does not depend on simplex_reduce anymore. I reimplemented simplex algorithm (see simplex_reduce_bounded function) to include upper bounding technique and faster pivoting. Now simple bounds such as l<=x<=u are not included in simplex tableau but handled separately, which is very beneficial when it comes to branch&bound as the tableau does not grow in size as bounding goes on. Also, basic columns are not kept in tableau anymore, only nonbasic are, which speeds up pivoting and saves up the space. I've adapted Bland's rule accordingly, so it is guaranteed that new algorithm does not cycle.
3. Branch&bound algorithm is now more controllable as it implements different (and more sophisticated) branching strategies. It is now capable to explore hundreds of thousands of nodes on moderate problem without hanging, slowing down or freezing. Inefficient tree generation using recursion is dropped, now only active nodes are stored in a list and tree exploration is done within while loop. lp_integertolerance was dropped in favor of better parameter doing the similar thing, lp_gaptolerance (relative integrality gap threshold).
4. I removed interior-point method implementation included in previous version. Now lpsolve uses GLPK solvers to operate in floating-point arithmetic. They are very fast. User can choose between primal simplex and interior-point algorithm.
5. Gomory cuts are now selected more efficiently and reoptimizations after appying cuts are much faster (see solve_relaxation function) as they start with the last feasible basis. Number of cuts that can be applied to a node can be limited or cuts can be disabled by using lp_maxcuts option.
6. Problems in CPLEX and MPS format can be loaded from files using GLPK library. All parameters and names are converted to Giac gen structures, so loading huge problems is disabled because they would eat up the resources. Loaded problems are treated as entered ones so user may apply additional bounds and integrality constraints to them before solving.

GLPK is a stable and maintained GNU library so it could be a new dependency for Giac without creating problems. But also I can use preprocessor directives to make its facilities optional before compilation. Bernard, do you think that it should be made optional? lpsolve can do much work on its own, however interior-point algorithm and loading from files would be disabled without GLPK.
Pièces jointes
lpsolve-fix5.zip
(25.05 Kio) Téléchargé 109 fois

parisse
Messages : 5082
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: linear programming

Thanks, I will have a look ASAP, probably tomorrow, but I can answer to your question about making GPLK optional: if it's not too much work, then please make it optional.

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

It's not a problem at all, GLPK facilities are treated as a separate module. I think I'll have time to do it today.

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

Here it is. To compile with GLPK one has to define HAVE_GLPK macro (somewhere during configuration, I guess).
If Giac is not compiled with GLPK, an error message will be printed if user tries to load a problem from file. When solving problem with approx coefficients, a warning message will be printed and builtin solver will be used after making all coefficients exact.
Pièces jointes
lpsolve-fix6.zip
(17.53 Kio) Téléchargé 101 fois

parisse
Messages : 5082
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: linear programming

Thank you! I'm updating configure.in to check for GLPK, as well as dispatch.h and lexer_tab_int.h with the new options.
What's INT_MAX ? I had to replace it with RAND_MAX to be able to compile, I guess it's the same.
I also wonder why you must #define _GLIBCXX_USE_CXX11_ABI 0 in lpsolve.h.

Code : Tout sélectionner

lpsolve(2x+y-z+4,[x<=1,y>=2,x+3y-z=2,2x-y+z<=8,-x+y<=5])
lpsolve(-4x-5y,[x+2y<=6,5x+4y<=20,0<=x,0<=y])
lpsolve(-7x+2y,[4x-12y<=20,-x+3y<=3],x=-5..5,y=0..inf,lp_maximize=true)
lpsolve(x-y-2z+3,[-3x-y+z<=3,2x-3y>=4z,x-z=y,x>=0,y<=0],lp_maximize)
lpsolve(-x-y,[y<=3x+1/2,y<=-5x+2],assume=lp_nonnegative)
lpsolve(x+y,[x<=8,-x+y<=4,-x+2y>=6,2x+y<=25,3x+y>=18,-x+2y>=6],assume=lp_nonnegative)
lpsolve(45.55x1+21.87x2,[1.64x1+2.67x2<=31.2,2.11x1+2.3x2>=13.9],assume=lp_nonnegative)
lpsolve(3x+4y,[x<=4,x+3y<=15,-x+2y>=5,x-y>=9,x+y=6],assume=lp_nonnegative,lp_maximize=true)
lpsolve(-6x+4y+z,[5x-10y<=20,2z-3y=6,-x+3y<=3],x=1..20,y=0..inf)
lpsolve(-5x-7y,[7x+y<=35,-x+3y<=6],assume=integer)
lpsolve(x+3y+3z,[x+3y+2z<=7,2x+2y+z<=11],assume=lp_nonnegative,lp_integervariables=[x,z],lp_maximize)
lpsolve(2x+5y,[3x-y=1,x-y<=5],assume=nonnegint)
lpsolve(x1+x2,[2x1+5x2<=16,6x1+5x2<=30],assume=nonnegint,lp_maximize)
lpsolve(8x1+11x2+6x3+4x4,[5x1+7x2+4x3+3x4<=14],assume=lp_binary,lp_maximize)

Should I update the list of examples?

lukamar
Messages : 294
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: linear programming

1. INT_MAX is defined in limits.h, it represents max value of signed int, it is at least 32767*(2^15-1). What is RAND_MAX? I suppose it's a very large integer, in which case it should do fine. INT_MAX represents integer "infinity", because INT_MAX of iterations/nodes/milliseconds will never be reached in reality.
2. I had to put #define _GLIBCXX_USE_CXX11_ABI 0 in lpsolve.h because without it gcc is giving me this error:

Code : Tout sélectionner

lpsolve.cc:1643: error: undefined reference to giac::gen2string[abi:cxx11](giac::gen const&)'
followed by

Code : Tout sélectionner

:-1: error: collect2: error: ld returned 1 exit status
That's happening in 11th line of _lpsolve function, where user entered string must be converted from gen to const char*. I think you may erase this instruction as you obviously don't have problems with gen2string function. Are we using different tools to compile? I'm using gcc and all standard libraries from Ubuntu 17.10 repositories, everything installed by default, nothing fancy. Instructions in Makefile to compile lpsolve.cc and test.cc (just main function used for testing) are:

Code : Tout sélectionner

g++ -Wall -g -pedantic -c lpsolve.cc
g++ -Wall -g -pedantic lpsolve.o test.cc -lgiac -lgmp -lglpk
I'm obviously not compatible with your settings as I experienced some strange segfaults which I told you about earlier (trying to use plus_inf etc.). Also, problems start to arise when trying to interchange string and gen (like this experience with gen2string). It looked to me like c++11 standard should be forced before compilation, but I tried that by including -std=c++11 option to the processor and it didn't work. For some reason defining _GLIBCXX_USE_CXX11_ABI solved the problem with gen2string.
3. All examples you listed should work fine. I suggest adding just one more to illustrate lp_verbose option:

Code : Tout sélectionner

lpsolve(x1+x2,[1867x1+1913x2=3618894],assume=nonnegint,lp_verbose=true)`