Le forum de XCAS

Xcas: un logiciel libre de calcul formel
Nous sommes actuellement le Mar Aoû 22, 2017 1:46 am

Heures au format UTC




Publier un nouveau sujet Répondre au sujet  [ 9 messages ] 
Auteur Message
 Sujet du message: linear programming
MessagePublié: Mar Aoû 08, 2017 8:35 pm 
Hors-ligne

Inscrit le: Ven Juin 30, 2017 9:55 am
Messages: 14
Localisation: Zagreb, Croatia
Hi Bernard,

I would like to contribute to the Giac development by providing the c++ code for function 'lpsolve' that emulates (almost fully) the function LPSolve from Maple's Optimization package. From the comment in misc.cc source file I see that an implementation of the two-phase simplex method is on the TODO list. I implemented it (upon already written 'simplex_reduce') alongside with the branch&bound method for solving (mixed) integer and binary LP problems. These are the main tools for 'lpsolve' (no active set or interior point methods were implemented).

With this addition (few hundreds lines of code), Giac would be able to solve great variety of practical problems, those which can be modeled as (integer) LP problems, Also, the compatibility with Maple would improve. As I am particulary interested in Giac having an usable LP solver, I decided to write one myself. I would be happy if you'd find time to review the code (it's GNU/GPL) and tell me if it qualifies for the testing version of Giac. Also, how do I submit the code? Attaching documents to messages in this forum seems to be disabled.


Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Mer Aoû 09, 2017 6:10 am 
Hors-ligne

Inscrit le: Mar Déc 20, 2005 4:02 pm
Messages: 3894
That would be really great!
I think you can attach zip files. Or enter the source in a message inside code delimiters, or email it to me!


Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Lun Aoû 14, 2017 5:16 pm 
Hors-ligne

Inscrit le: Ven Juin 30, 2017 9:55 am
Messages: 14
Localisation: Zagreb, Croatia
Here is the code. Sorry for the delay, I had to fix some bugs.
The enclosed Makefile creates executable a.out which can be run to test the function lpsolve. Just enter the problem as "(obj,constr,....)". There are several examples of LP problems included, mostly from the Maplesoft web page explaining their function LPSolve. They should all work except the last, eleventh example in the documentation entry for function lpsolve, I have no idea why it does not work.
There are couple of things I should mention regarding lpsolve.cc. As the inclusion of 'plus_inf' and 'minus_inf' constants in the code made it segfault for some reason, I used a workaround by adding function 'make_infinity' which creates +infinity symbol. It should be removed and all instances of 'inifnity' outside the comments replaced by 'plus_inf', as in other Giac functions. Also, there was a bug in function 'simplex_reduce' which I had to correct, hence the fixed version 'simplex_reduce2'. When determining basis after pivoting, the 'simplex_reduce' only checked for zeros in an idn column, but not that the nonzero element is equal to one. I added that check (see line 171 in lpsolve.cc), which should be included in the code for 'simplex_reduce'.
For example, to solve the LP problem included in cascmd_fr (when explaining 'simplex_reduce' function), one should run a.out and enter:
Code:
(2x+y-z+4,[x<=1,y>=2,x+3y-z=2,2x-y+z<=8,-x+y<=5])

Also, the function 'lpsolve' requires a few new option keywords to be added to Giac. These are listed in lpsolve.h. I used the available option 'integer', but other options, such as 'nonnegative' and so on, should be defined in similar fashion. Since they're missing, use the corresponding integer values defined in lpsolve.h. For example, to specify 'assume=nonnegative' one should type '0=7'.


Pièces jointes:
lpsolve.zip [15.56 Kio]
Téléchargé 2 fois
Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Mar Aoû 15, 2017 1:58 pm 
Hors-ligne

Inscrit le: Mar Déc 20, 2005 4:02 pm
Messages: 3894
I have mirrored your bugfix in misc.cc, so that you can remove simplex_reduce2.
For the lp keywords, I can add them to the lexer, except that:
* assume is already a giac keyword, you can't parse it as an integer, but it's not difficult to check for it.
* nonnegint. already exists, value is _NONNEGINT=4*256+2,
I can add the lp specific integer constants to maple_enum in dispatch.h, values may therefore change with respect to lpsolve.h (lpsolve.cc unchanged).
This could bring something like:
Code:
      if (R==at_equal) {
        // parse the argument in form "option=value"
        vecteur ops(*((*it)._SYMBptr->feuille._VECTptr));
        if (ops[0].type==_IDNT && intrv2vec(ops[1],lr,contextptr)) {
          // the boundaries for variable ops[0] are set
          if (!is_zero(lr[0]) && is_strictly_greater(lr[0],-infinity,contextptr))
            bconstr.push_back(symbolic(at_superieur_egal,makevecteur(ops[0],lr[0])));
          if (is_strictly_greater(infinity,lr[1],contextptr))
            bconstr.push_back(symbolic(at_inferieur_egal,makevecteur(ops[0],lr[1])));
          if (is_positive(lr[0],contextptr))
            nvars.push_back(ops[0]);
          continue;
        }
        if (ops[0]==at_assume){ // put here the case LP_ASSUME:
          ...
         continue;
        }
        if (is_integer(ops[0])){
          ...
        }
     }

Once all is fixed, I can simply add lpsolve.cc/lpsolve.h to Makefile.am and lpsolve should be available in Giac/Xcas. I don't know why you get segfault with the builtin inf, perhaps a loader order bug, anyway this should disappear once the code is integrated in the library.


Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Mar Aoû 15, 2017 2:21 pm 
Hors-ligne

Inscrit le: Mar Déc 20, 2005 4:02 pm
Messages: 3894
Actually, I think it's safer to add the prefix lp_ to all constants that the user will key in (while accepting both lp_integer and integer or lp_nonnegint and nonnegint should be supported). Code added in lexer_tab_int.h:
Code:
      {"lp_binary"               ,0, _LP_BINARY, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_binaryvariables"               ,0, _LP_BINARYVARIABLES, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_depthlimit"               ,0, _LP_DEPTHLIMIT, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_integer"               ,0, _LP_INTEGER, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_integervariables"               ,0, _LP_INTEGERVARIABLES, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_maximize"               ,0, _LP_MAXIMIZE, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_nonnegative"               ,0, _LP_NONNEGATIVE, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_nonnegint"               ,0, _LP_NONNEGINT, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_variables"               ,0, _LP_VARIABLES, _INT_MAPLECONVERSION, T_TYPE_ID},

In dispatch.h:
Code:
  enum maple_conversion {
    _TRIG=100,
    _EXPLN=101,
    _PARFRAC=102,
    _FULLPARFRAC=103,
    _BASE=104,
    _CONFRAC=105,
    _MAPLE_LIST=256,
    _POSINT=1*256+2,
    _NEGINT=2*256+2,
    _NONPOSINT=3*256+2,
    _NONNEGINT=4*256+2,
    _LP_BINARY=106,           //   * binary
    _LP_BINARYVARIABLES=107,  //   * binaryvariables
    _LP_DEPTHLIMIT=108,       //   * depthlimit
    _LP_INTEGER=109,          //   * integer
    _LP_INTEGERVARIABLES=110, //   * integervariables
    _LP_MAXIMIZE=111,         //   * maximize
    _LP_NONNEGATIVE=112,      //   * nonnegative
    _LP_NONNEGINT=113,        //   * nonnegint
    _LP_VARIABLES=114,        //   * variables
  };

Code also added in gen.cc in printint32.
Code:
    ...
      case _NONNEGINT:
   return "nonnegint";
      case _LP_BINARY:
   return "lp_binary";
      case _LP_BINARYVARIABLES:
   return "lp_binaryvariables";
      case _LP_DEPTHLIMIT:
   return "lp_depthlimit";
      case _LP_MAXIMIZE:
   return "lp_maximize";
      case _LP_NONNEGATIVE:
   return "lp_nonnegative";
      case _LP_NONNEGINT:
   return "lp_nonnegint";
      case _LP_VARIABLES:
   return "lp_variables";
      case _LP_INTEGER:
   return "lp_integer";
      case _LP_INTEGERVARIABLES:
   return "lp_integervariables";
      }


Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Mer Aoû 16, 2017 2:32 pm 
Hors-ligne

Inscrit le: Ven Juin 30, 2017 9:55 am
Messages: 14
Localisation: Zagreb, Croatia
Here is the final version of the code. I fixed the bug which caused wrong result to the example 11, it was in branch_and_bound function when comparing with the incumbent solution. Now all examples should work. I put plus_inf and minus_inf on the respective places and the code compiles fine. Also, i made a few improvements, did some refactoring, added more examples and updated the documentation blocks.
I also plan to write an entry to cascmd TeX file (english version) about the function 'lpsolve'. I'll send it to you when I finish, most probably in the next day or two.

Edit: now I see that I forgot to change the names in lpsolve.h. lpmax should be simplex_twophase, lpmin should be simplex_dual and ilpbb should be branch_and_bound.


Pièces jointes:
lpsolve-fixed.zip [13.83 Kio]
Téléchargé 2 fois
Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Dim Aoû 20, 2017 5:57 pm 
Hors-ligne

Inscrit le: Mar Déc 20, 2005 4:02 pm
Messages: 3894
It is now included in the javascript version if you want to check:
https://www-fourier.ujf-grenoble.fr/~parisse/xcasfr.html
I will update giac source and version during the first week of September (perhaps before if I can).


Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Dim Aoû 20, 2017 7:32 pm 
Hors-ligne

Inscrit le: Ven Juin 30, 2017 9:55 am
Messages: 14
Localisation: Zagreb, Croatia
Thanks. But the web interface does not recognize the 'lpsolve' command. I cleared the cached web content and reloaded, but with no effect. How do I update CAS? I'm using Firefox.


Haut
 Profil  
 
 Sujet du message: Re: linear programming
MessagePublié: Lun Aoû 21, 2017 5:51 am 
Hors-ligne

Inscrit le: Mar Déc 20, 2005 4:02 pm
Messages: 3894
The UI does not display lpsolve as a command (I must add it in a file somewhere) but it should work.


Haut
 Profil  
 
Afficher les messages publiés depuis:  Trier par  
Publier un nouveau sujet Répondre au sujet  [ 9 messages ] 

Heures au format UTC


Qui est en ligne ?

Utilisateurs parcourant actuellement ce forum : Aucun utilisateur inscrit et 1 invité


Vous ne pouvez pas publier de nouveaux sujets dans ce forum
Vous ne pouvez pas répondre aux sujets dans ce forum
Vous ne pouvez pas éditer vos messages dans ce forum
Vous ne pouvez pas supprimer vos messages dans ce forum
Vous ne pouvez pas insérer de pièces jointes dans ce forum

Rechercher pour:
Sauter vers:  
cron
Powered by phpBB® Forum Software © phpBB Group
Traduction réalisée par Maël Soucaze © 2009 phpBB.fr