1 /*
2 * lpsolve.cc
3 *
4 * Copyright 2017 Luka Marohnić
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 #include "giacPCH.h"
21 #include "giac.h"
22 #include "lpsolve.h"
23 #include "optimization.h"
24 #include <ctime>
25
26 #ifndef DBL_MAX
27 #define DBL_MAX 1.79769313486e+308
28 #endif
29
30 using namespace std;
31
32 #ifndef NO_NAMESPACE_GIAC
33 namespace giac {
34 #endif //ndef NO_NAMESPACE_GIAC
35
36 const gen inf=symbolic(at_plus,unsigned_inf);
37 const gen minf=symbolic(at_neg,unsigned_inf);
38
39 /*
40 * Convert double to exact gen.
41 */
double2fraction(double d,GIAC_CONTEXT)42 gen double2fraction(double d,GIAC_CONTEXT) {
43 return exact(double2gen(d),contextptr);
44 }
45
46 /*
47 * Convert gen to double.
48 */
gen2double(const gen & g,GIAC_CONTEXT)49 double gen2double(const gen &g,GIAC_CONTEXT) {
50 return _evalf(g,contextptr).DOUBLE_val();
51 }
52
53 /*
54 * Return number of digits of an unsigned integer.
55 */
numdigits(unsigned i)56 int numdigits(unsigned i) {
57 return 1+(i>0?(int)std::log10((double)i):0);
58 }
59
60 /*
61 * Return true iff g is a (vector of) real number or +infinity or -infinity.
62 */
is_realcons(const gen & g,GIAC_CONTEXT)63 bool is_realcons(const gen &g,GIAC_CONTEXT) {
64 if (g.type==_VECT) {
65 vecteur &v = *g._VECTptr;
66 for (const_iterateur it=v.begin();it!=v.end();++it) {
67 if (!is_realcons(*it,contextptr))
68 return false;
69 }
70 return true;
71 }
72 return (is_inf(g) || _evalf(g,contextptr).type==_DOUBLE_);
73 }
74
75 /*
76 * If g is an interval, store its bounds to pair p and return true. If g is not
77 * an interval, return false.
78 */
interval2pair(const gen & g,pair<gen,gen> & p,GIAC_CONTEXT)79 bool interval2pair(const gen &g,pair<gen,gen> &p,GIAC_CONTEXT) {
80 if (g.type!=_SYMB || !g.is_symb_of_sommet(at_interval))
81 return false; //g is not an interval
82 vecteur &v=*g._SYMBptr->feuille._VECTptr;
83 p=make_pair(v[0],v[1]);
84 return is_realcons(v,contextptr);
85 }
86
87 /*
88 * Make singleton vector (with one at position j and zeros at other positions).
89 */
singleton(int n,int j)90 vecteur singleton(int n,int j) {
91 vecteur v(n,gen(0));
92 v[j]=gen(1);
93 return v;
94 }
95
96 /*
97 * Insert column in matrix at position j.
98 */
insert_column(matrice & m,const vecteur & c,int j)99 void insert_column(matrice &m,const vecteur &c,int j) {
100 assert(m.size()==c.size());
101 for (int i=0;i<int(m.size());++i) {
102 m[i]._VECTptr->insert(j>=0?m[i]._VECTptr->begin()+j:m[i]._VECTptr->end()+j,c[i]);
103 }
104 }
105
106 /*
107 * Append column to matrix.
108 */
append_column(matrice & m,const vecteur & c)109 void append_column(matrice &m,const vecteur &c) {
110 assert(m.size()==c.size());
111 matrice mt=mtran(m);
112 mt.push_back(c);
113 m=mtran(mt);
114 }
115
116 /*
117 * Remove jth column from matrix. If j<0, count from last column towards the first.
118 */
remove_column(matrice & m,int j)119 void remove_column(matrice &m,int j) {
120 for (int i=0;i<int(m.size());++i) {
121 m[i]._VECTptr->erase(j+(j>=0?m[i]._VECTptr->begin():m.end()));
122 }
123 }
124
125 /*
126 * Get jth column from matrix. If j<0, count from last column towards the first.
127 */
jth_column(const matrice & m,int j)128 vecteur jth_column(const matrice &m,int j) {
129 int n=m.front()._VECTptr->size();
130 vecteur col(m.size());
131 for (int i=0;i<int(m.size());++i) {
132 col[i]=m[i][j>=0?j:n+j];
133 }
134 return col;
135 }
136
137 /*
138 * Multiply coefficients by LCM of denominators and then divide them by their
139 * GCD.
140 */
integralize(const vecteur & v_orig,GIAC_CONTEXT)141 vecteur integralize(const vecteur &v_orig,GIAC_CONTEXT) {
142 vecteur v(v_orig),vd;
143 for (const_iterateur it=v.begin();it!=v.end();++it) {
144 if (!is_zero(*it))
145 vd.push_back(_denom(*it,contextptr));
146 }
147 if (vd.empty())
148 return v;
149 v=multvecteur(abs(_lcm(vd,contextptr),contextptr),v);
150 return divvecteur(v,abs(_gcd(v,contextptr),contextptr));
151 }
152
153 /*
154 * LP variable constructor. By default, it is a nonnegative variable
155 * unrestricted from above.
156 */
lp_variable()157 lp_variable::lp_variable() {
158 is_integral=false;
159 sign_type=_LP_VARSIGN_POS;
160 range=lp_range();
161 range.lbound=gen(0);
162 fill_n(nbranch,2,0);
163 }
164
165 /*
166 * Update lower (dir=0) or upper (dir=1) pseudocost.
167 */
update_pseudocost(double delta,double fr,int dir)168 void lp_variable::update_pseudocost(double delta, double fr, int dir) {
169 double sigma=pseudocost[dir]*nbranch[dir];
170 sigma+=delta/(dir==0?fr:(1-fr));
171 pseudocost[dir]=sigma/(++nbranch[dir]);
172 }
173
174 /*
175 * Return score, a positive value based on pseudocost values. Variable with the
176 * best (highest) score is selected for branching.
177 */
score(double fr)178 double lp_variable::score(double fr) {
179 if (nbranch[0]==0 || nbranch[1]==0)
180 return 0;
181 double qlo=fr*pseudocost[0],qhi=(1-fr)*pseudocost[1];
182 return (1-LP_SCORE_FACTOR)*std::min(qlo,qhi)+LP_SCORE_FACTOR*std::max(qlo,qhi);
183 }
184
185 /*
186 * Range constructor: by default, it contains no restriction.
187 */
lp_range()188 lp_range::lp_range () {
189 lbound=minf;
190 ubound=inf;
191 }
192
193 /*
194 * Settings constructor loads some sensible defaults.
195 */
lp_settings()196 lp_settings::lp_settings() {
197 verbose=false;
198 status_report_freq=0.2;
199 solver=_LP_SIMPLEX;
200 precision=_LP_PROB_DEPENDENT;
201 presolve=false;
202 maximize=false;
203 relative_gap_tolerance=0.0;
204 has_binary_vars=false;
205 varselect=-1;
206 nodeselect=_LP_BEST_PROJECTION;
207 depth_limit=RAND_MAX;
208 node_limit=RAND_MAX;
209 iteration_limit=RAND_MAX;
210 time_limit=RAND_MAX;
211 max_cuts=5;
212 }
213
214 /*
215 * Stats constructor initializes the status container for the problem being
216 * solved. It is used to monitor the progress and to summarize when done.
217 */
lp_stats()218 lp_stats::lp_stats() {
219 subproblems_examined=0;
220 cuts_applied=0;
221 max_active_nodes=0;
222 mip_gap=-1; //negative means undefined
223 }
224
225 /*
226 * Pivot on element with coordinates I,J.
227 */
pivot_ij(matrice & m,int I,int J,bool negate=false)228 void pivot_ij(matrice &m,int I,int J,bool negate=false) {
229 int nr=m.size(),nc=m.front()._VECTptr->size();
230 gen a(m[I][J]),b;
231 m[I]=divvecteur(*m[I]._VECTptr,a);
232 vecteur col(jth_column(m,J)),&pv=*m[I]._VECTptr;
233 for (int i=0;i<nr;++i) {
234 m[i]._VECTptr->at(J)=(i==I?gen(negate?-1:1)/a:gen(0));
235 }
236 for (int i=0;i<nr;++i) {
237 b=col[i];
238 if (i==I || is_zero(b))
239 continue;
240 vecteur &v=*m[i]._VECTptr;
241 for (int j=0;j<nc;++j) {
242 v[j]-=b*pv[j];
243 }
244 }
245 }
246
247 /*
248 * Simplex algorithm that handles upper bounds of the variables. The solution x
249 * satisfies 0<=x<=u. An initial basis must be provided.
250 *
251 * Basis is a vector of integers B and B[i]=j means that jth variable is basic
252 * and appears in ith row (constraint). Basic columns are not kept in matrix,
253 * which contains only the columns of nonbasic variables. A nonbasic variable
254 * is assigned to the respective column with integer vector 'cols': cols[i]=j
255 * means that ith column of the matrix is associated with the jth variable. The
256 * algorithm uses upper-bounding technique when pivoting and uses an adaptation
257 * of Bland's rule to prevent cycling. ith element of 'is_slack' is true iff
258 * ith (nonbasic) variable is at its upper bound.
259 *
260 * If limit>0, simplex algorithm will terminate after that many iterations.
261 */
simplex_reduce_bounded(const matrice & m_orig,const vecteur & u,vector<bool> & is_slack,vecteur & bfs,gen & optimum,ints & basis,ints & cols,int limit,GIAC_CONTEXT)262 matrice simplex_reduce_bounded(const matrice &m_orig,const vecteur &u,vector<bool> &is_slack,vecteur &bfs,gen &optimum,
263 ints &basis,ints &cols,int limit,GIAC_CONTEXT) {
264 matrice m(m_orig);
265 //ev, lv: indices of entering and leaving variables
266 //ec, lr: 'entering' column and 'leaving' row in matrix m, respectively
267 int nr=basis.size(),nc=cols.size(),nv=nr+nc,ec,ev,lr,lv,icount=0;
268 gen a,b,ratio,mincoeff,opt;
269 //iterate the simplex method
270 bool choose_first=false;
271 optimum=undef;
272 while ((icount++)<limit) {
273 opt=m.back()._VECTptr->back();
274 //determine which variable enters the basis
275 mincoeff=gen(0);
276 ev=-1;
277 m.back()=_eval(m.back(),contextptr);
278 vecteur &last=*m.back()._VECTptr;
279 for (int j=0;j<nc;++j) {
280 int k=cols[j];
281 if ((choose_first && is_strictly_positive(-last[j],contextptr) &&
282 (ev<0 || k+(is_slack[k]?nv:0)<ev+(is_slack[ev]?nv:0))) ||
283 (!choose_first && is_strictly_greater(mincoeff,last[j],contextptr))) {
284 ec=j;
285 ev=k;
286 mincoeff=last[j];
287 }
288 }
289 if (ev<0) //the current solution is optimal
290 break;
291 //determine which variable leaves the basis
292 mincoeff=inf;
293 lv=-1;
294 bool hits_ub,ub_subs;
295 for (int i=0;i<nr;++i) {
296 a=m[i][ec];
297 b=m[i]._VECTptr->back();
298 int j=basis[i];
299 if (is_strictly_positive(a,contextptr) && is_greater(mincoeff,ratio=b/a,contextptr))
300 hits_ub=false;
301 else if (is_strictly_positive(-a,contextptr) && !is_inf(u[j]) &&
302 is_greater(mincoeff,ratio=(b-u[j])/a,contextptr))
303 hits_ub=true;
304 else continue;
305 if (is_strictly_greater(mincoeff,ratio,contextptr)) {
306 lv=-1;
307 mincoeff=ratio;
308 }
309 if (lv<0 || (choose_first && j+(is_slack[j]?nv:0)<lv+(is_slack[lv]?nv:0))) {
310 lv=j;
311 lr=i;
312 ub_subs=hits_ub;
313 }
314 }
315 if (lv<0 && is_inf(u[ev])) { //solution is unbounded
316 optimum=inf;
317 return m;
318 }
319 if (is_zero(mincoeff))
320 choose_first=true; //switch on Bland's rule
321 if (lv<0 || is_greater(mincoeff,u[ev],contextptr)) {
322 for (iterateur it=m.begin();it!=m.end();++it) {
323 it->_VECTptr->back()-=u[ev]*(a=it->_VECTptr->at(ec));
324 it->_VECTptr->at(ec)=-a;
325 }
326 if (ev<nv)
327 is_slack[ev]=!is_slack[ev];
328 continue;
329 }
330 //swap variables: basic leaves, nonbasic enters
331 if (ub_subs) {
332 m[lr]._VECTptr->back()-=u[lv];
333 if (lv<nv)
334 is_slack[lv]=!is_slack[lv];
335 }
336 pivot_ij(m,lr,ec,ub_subs);
337 basis[lr]=ev;
338 cols[ec]=lv;
339 }
340 m=*_eval(m,contextptr)._VECTptr;
341 optimum=m[nr][nc];
342 bfs=vecteur(nv,gen(0));
343 for (int i=0;i<nr;++i) {
344 bfs[basis[i]]=m[i][nc];
345 }
346 for (int j=0;j<nv;++j) {
347 if (is_slack[j])
348 bfs[j]=u[j]-bfs[j];
349 }
350 return m;
351 }
352
353 /*
354 * Solve the relaxed subproblem corresponding to this node.
355 *
356 * This function uses two-phase simplex method and applies suitable Gomory
357 * mixed integer cuts generated after (each re)optimization. Weak GMI cuts are
358 * discarded either because of small away or because not being parallel enough
359 * to the objective. Cuts with too large coefficients (when integralized) are
360 * discarded too because they slow down the computational process. Generated
361 * cuts are kept in the problem structure to be used by child suboproblems
362 * during the branch&bound algorithm.
363 */
solve_relaxation()364 int lp_node::solve_relaxation() {
365 int nrows=prob->constr.nrows(),ncols=prob->constr.ncols(),bs;
366 matrice m;
367 vecteur obj(prob->objective.first),l(ncols),u(ncols),br,row,b,lh,gmi_cut;
368 gen rh,obj_ct(0),mgn,minmgn;
369 ints cols(ncols),basis;
370 vector<bool> is_slack(ncols,false);
371 map<int,int> slack_cut;
372 bool is_mip=prob->has_integral_variables();
373 //determine upper and lower bounds
374 for (int j=0;j<ncols;++j) {
375 lp_variable &var=prob->variables[j];
376 lp_range &range=ranges[j];
377 l[j]=max(var.range.lbound,range.lbound,prob->ctx);
378 u[j]=min(var.range.ubound,range.ubound,prob->ctx);
379 if (is_strictly_greater(l[j],u[j],prob->ctx))
380 return _LP_INFEASIBLE;
381 }
382 //populate matrix with consraint coefficients
383 m=*_matrix(makesequence(nrows,ncols+1,0),prob->ctx)._VECTptr;
384 for (int i=0;i<nrows;++i) for (int j=0;j<ncols;++j) {
385 m[i]._VECTptr->at(j)=prob->constr.lhs[i][j];
386 }
387 b=prob->constr.rhs;
388 //shift variables according to their lower bounds such that l<=x<=u becomes 0<=x'<=u'
389 for (int j=0;j<ncols;++j) {
390 b=subvecteur(b,multvecteur(l[j],jth_column(m,j)));
391 u[j]-=l[j];
392 obj_ct+=obj[j]*l[j];
393 cols[j]=j;
394 }
395 for (int i=0;i<nrows;++i) {
396 m[i]._VECTptr->at(ncols)=b[i];
397 }
398 //assure that the right-hand side column has nonnegative coefficients
399 for (iterateur it=m.begin();it!=m.end();++it) {
400 if (is_strictly_positive(-it->_VECTptr->back(),prob->ctx)) {
401 *it=-*it;
402 }
403 }
404 //append cuts inherited from parent node
405 for (ints::const_iterator it=cut_indices.begin();it!=cut_indices.end();++it) {
406 insert_column(m,vecteur(nrows,gen(0)),-1);
407 slack_cut[ncols]=*it;
408 prob->cuts.get_lr(*it,lh,rh);
409 for (int i=0;i<prob->nv();++i) {
410 rh-=lh[i]*l[i];
411 }
412 row=vecteur(cols.size(),gen(0));
413 for (int i=0;i<int(cols.size());++i) {
414 int j=cols.at(i);
415 if (j<prob->nv())
416 row[i]=lh[j];
417 }
418 row.push_back(gen(-1));
419 row.push_back(rh);
420 if (is_strictly_positive(-rh,prob->ctx))
421 row=multvecteur(gen(-1),row);
422 m.push_back(row);
423 l.push_back(gen(0));
424 u.push_back(inf);
425 is_slack.push_back(false);
426 obj.push_back(gen(0));
427 cols.push_back(ncols++);
428 ++nrows;
429 }
430 //optimize-add cut-reoptimize-add cut...
431 //repeat until no more cuts are generated or max_cuts limit is reached
432 while (true) {
433 br=vecteur(int(cols.size())+1,gen(0));
434 bs=basis.size();
435 basis.resize(nrows);
436 u.resize(ncols+nrows-bs);
437 for (int i=bs;i<nrows;++i) {
438 br=subvecteur(br,*m[i]._VECTptr);
439 basis[i]=ncols+i;
440 u[ncols+i-bs]=inf;
441 }
442 m.push_back(br);
443 //phase 1: minimize the sum of artificial variables
444 m=simplex_reduce_bounded(m,u,is_slack,solution,optimum,basis,cols,RAND_MAX,prob->ctx);
445 if (!is_zero(_simplify(optimum,prob->ctx)))
446 return _LP_INFEASIBLE; //at least one artificial variable is basic and positive
447 m.pop_back(); //remove bottom row
448 for (int i=0;i<nrows;++i) {
449 int j=basis[i];
450 if (j<ncols)
451 continue;
452 //ith basic variable is artificial, push it out of the basis
453 int k=0;
454 for (;k<ncols && (is_zero(m[i][k]) || cols[k]>=ncols);++k);
455 if (k==ncols)
456 return _LP_ERROR;
457 pivot_ij(m,i,k);
458 basis[i]=cols[k];
459 cols[k]=j;
460 }
461 //phase 1 finished, remove artificial columns from m
462 for (int j=int(cols.size())-1;j>=0;--j) {
463 if (cols[j]>=ncols) {
464 remove_column(m,j);
465 cols.erase(cols.begin()+j);
466 }
467 }
468 //append bottom row to maximize -obj
469 br=vecteur(ncols-nrows+1);
470 for (int j=0;j<ncols-nrows;++j) {
471 int k=cols[j];
472 br[j]=obj[k];
473 if (is_slack[k]) {
474 br.back()-=br[j]*u[k];
475 br[j]=-br[j];
476 }
477 }
478 for (int i=0;i<nrows;++i) {
479 int j=basis[i];
480 if (is_slack[j])
481 br.back()-=obj[j]*u[j];
482 br=subvecteur(br,multvecteur(is_slack[j]?-obj[j]:obj[j],*m[i]._VECTptr));
483 }
484 m.push_back(br);
485 u.resize(ncols);
486 //phase 2: optimize the objective
487 m=simplex_reduce_bounded(m,u,is_slack,solution,optimum,basis,cols,
488 prob->settings.iteration_limit,prob->ctx);
489 if (is_inf(optimum))
490 return _LP_UNBOUNDED; //solution is unbounded
491 m.pop_back(); //remove bottom row
492 if (!is_mip) break;
493 if (int(cut_indices.size())>=prob->settings.max_cuts)
494 break;
495 //try to generate Gomory mixed integer cut
496 gmi_cut.clear();
497 for (int i=0;i<nrows;++i) {
498 gen p(fracpart(solution[basis[i]]));
499 vecteur eq(*m[i]._VECTptr);
500 if (!is_zero(p) && !is_integer(eq.back())) {
501 gen f0(fracpart(eq.back())),fj,sp(0),eqnorm(0);
502 double away=_evalf(min(f0,gen(1)-f0,prob->ctx),prob->ctx).DOUBLE_val();
503 if (away<LP_MIN_AWAY)
504 continue; //too small away, discard this cut
505 eq.back()=gen(1);
506 for (int k=0;k<int(cols.size());++k) {
507 int j=cols[k];
508 if (j<prob->nv() && prob->variables[j].is_integral) {
509 fj=fracpart(eq[k]);
510 eq[k]=(is_strictly_greater(fj,f0,prob->ctx)?(fj-1)/(f0-1):fj/f0);
511 }
512 else
513 eq[k]=eq[k]/(is_strictly_positive(eq[k],prob->ctx)?f0:(f0-1));
514 if (j<prob->nv()) {
515 sp+=eq[k]*obj[j];
516 eqnorm+=pow(eq[k],2);
517 }
518 }
519 double dsp=_evalf(sp,prob->ctx).DOUBLE_val(),deqnorm=_evalf(eqnorm,prob->ctx).DOUBLE_val();
520 if (std::abs(dsp/(prob->objective_norm*std::sqrt(deqnorm)))<LP_MIN_PARALLELISM)
521 continue; //this cut is not parallel enough to the objective
522 eq=integralize(eq,prob->ctx); //turn cut into an equivalent integral representation
523 mgn=_max(_abs(eq,prob->ctx),prob->ctx);
524 if (_evalf(mgn,prob->ctx).DOUBLE_val()>LP_MAX_MAGNITUDE)
525 continue; //this cut has too large coefficients
526 if (gmi_cut.empty() || is_strictly_greater(minmgn,mgn,prob->ctx)) {
527 minmgn=mgn;
528 gmi_cut=eq;
529 }
530 }
531 }
532 if (gmi_cut.empty())
533 break; //no acceptable cut was generated, so there's no need to reoptimize further
534 //store GMI cut
535 lh=vecteur(ncols,gen(0));
536 rh=gmi_cut.back();
537 for (int k=0;k<int(cols.size());++k) {
538 int j=cols[k];
539 if(!is_zero(lh[j]=gmi_cut[k])) {
540 if (is_slack[j]) {
541 rh-=u[j]*lh[j];
542 lh[j]=-lh[j];
543 }
544 rh+=l[j]*lh[j];
545 }
546 }
547 vecteur slacks(lh.begin()+prob->nv(),lh.end()),orig_lh;
548 gen orig_rh;
549 lh.resize(prob->nv());
550 for (int k=0;k<int(slacks.size());++k) {
551 int j=prob->nv()+k;
552 prob->cuts.get_lr(slack_cut[j],orig_lh,orig_rh);
553 rh+=slacks[k]*orig_rh;
554 lh=addvecteur(lh,multvecteur(slacks[k],orig_lh));
555 }
556 cut_indices.push_back(prob->cuts.nrows());
557 prob->cuts.append(lh,rh,_LP_GEQ);
558 ++prob->stats.cuts_applied;
559 //append GMI cut to simplex tableau and reoptimize
560 gmi_cut.insert(gmi_cut.end()-1,gen(-1));
561 insert_column(m,vecteur(nrows,gen(0)),-1);
562 m.push_back(gmi_cut);
563 slack_cut[ncols]=cut_indices.back();
564 l.push_back(gen(0));
565 u.push_back(inf);
566 is_slack.push_back(false);
567 obj.push_back(gen(0));
568 cols.push_back(ncols++);
569 ++nrows;
570 }
571 for (int i=0;i<prob->nv();++i) {
572 solution[i]+=l[i];
573 }
574 solution.resize(prob->nv());
575 optimum=obj_ct-optimum;
576 //compute some data useful during branch&bound
577 opt_approx=gen2double(optimum,prob->ctx);
578 infeas=gen(0);
579 most_fractional=-1;
580 gen p,ifs,max_ifs(0);
581 for (int i=0;i<prob->nv();++i) {
582 if (!prob->variables[i].is_integral)
583 continue;
584 p=fracpart(solution[i]);
585 ifs=min(p,gen(1)-p,prob->ctx);
586 if (is_zero(ifs))
587 continue;
588 fractional_vars[i]=gen2double(p,prob->ctx);
589 infeas+=ifs;
590 if (is_strictly_greater(ifs,max_ifs,prob->ctx)) {
591 most_fractional=i;
592 max_ifs=ifs;
593 }
594 }
595 return _LP_SOLVED;
596 }
597
598 /*
599 * Return fractional part of g, i.e. [g]=g-floor(g). It is always 0<=[g]<1.
600 */
fracpart(const gen & g)601 gen lp_node::fracpart(const gen &g) {
602 return g-_floor(g,prob->ctx);
603 }
604
605 /*
606 * Return true iff variable with specified index is fractional.
607 */
is_var_fractional(int index)608 bool lp_node::is_var_fractional(int index) {
609 for (map<int,double>::const_iterator it=fractional_vars.begin();it!=fractional_vars.end();++it) {
610 if (it->first==index)
611 return true;
612 }
613 return false;
614 }
615
616 /*
617 * Create a child node with copy of ranges and cut indices and depth increased
618 * by one.
619 */
create_child()620 lp_node lp_node::create_child() {
621 lp_node node;
622 node.depth=depth+1;
623 node.ranges=vector<lp_range>(ranges);
624 node.cut_indices=ints(cut_indices);
625 node.prob=prob;
626 return node;
627 }
628
629 /*
630 * Return true iff there are integrality constraints, i.e. iff this is a
631 * (mixed) integer problem.
632 */
has_integral_variables()633 bool lp_problem::has_integral_variables() {
634 for (vector<lp_variable>::const_iterator it=variables.begin();it!=variables.end();++it) {
635 if (it->is_integral)
636 return true;
637 }
638 return false;
639 }
640
641 /*
642 * Return true iff problem has approximate (floating-point) coefficients.
643 */
has_approx_coefficients()644 bool lp_problem::has_approx_coefficients() {
645 if (is_approx(objective.first) ||
646 objective.second.is_approx() ||
647 is_approx(constr.lhs) ||
648 is_approx(constr.rhs))
649 return true;
650 for (vector<lp_variable>::const_iterator it=variables.begin();it!=variables.end();++it) {
651 if (it->range.lbound.is_approx() || it->range.ubound.is_approx())
652 return true;
653 }
654 return false;
655 }
656
657 /*
658 * Set objective function parameters.
659 */
set_objective(const vecteur & v,const gen & ft)660 void lp_problem::set_objective(const vecteur &v, const gen &ft) {
661 objective.first=v;
662 objective.second=ft;
663 for (const_iterateur it=v.begin();it!=v.end();++it) {
664 obj_approx.push_back(gen2double(abs(*it,ctx),ctx));
665 }
666 }
667
668 /*
669 * Display a message.
670 */
message(const char * msg,bool err)671 void lp_problem::message(const char *msg, bool err) {
672 if (err || settings.verbose)
673 *logptr(ctx) << msg << "\n";
674 }
675
676 /*
677 * Return the jth column (coefficients of the jth variable) as vecteur.
678 */
column(int index)679 vecteur lp_constraints::column(int index) {
680 return jth_column(lhs,index);
681 }
682
683 /*
684 * Duplicate the jth column.
685 */
duplicate_column(int index)686 void lp_constraints::duplicate_column(int index) {
687 assert(index<ncols());
688 vecteur col(column(index));
689 insert_column(lhs,col,index);
690 }
691
692 /*
693 * Multiply the jth column by -1.
694 */
negate_column(int index)695 void lp_constraints::negate_column(int index) {
696 for (int i=0;i<nrows();++i) {
697 vecteur &lh=*lhs[i]._VECTptr;
698 lh[index]=-lh[index];
699 }
700 }
701
702 /*
703 * Subtract v from rhs column of constraints.
704 */
subtract_from_rhs_column(const vecteur & v)705 void lp_constraints::subtract_from_rhs_column(const vecteur &v) {
706 assert(int(v.size())==nrows());
707 for (int i=0;i<nrows();++i) {
708 rhs[i]-=v[i];
709 }
710 }
711
712 /*
713 * Append constraint "lh rel rh".
714 */
append(const vecteur & lh,const gen & rh,int relation_type)715 void lp_constraints::append(const vecteur &lh, const gen &rh, int relation_type) {
716 assert(nrows()==0 || int(lh.size())==ncols());
717 lhs.push_back(lh);
718 rhs.push_back(rh);
719 rv.push_back(relation_type);
720 }
721
722 /*
723 * Set the constraint with specified index.
724 */
set(int index,const vecteur & lh,const gen & rh,int relation_type)725 void lp_constraints::set(int index, const vecteur &lh, const gen &rh, int relation_type) {
726 assert(index<nrows());
727 lhs[index]=lh;
728 rhs[index]=rh;
729 rv[index]=relation_type;
730 }
731
732 /*
733 * Get left and right side of the constraint with specified index.
734 */
get_lr(int index,vecteur & lh,gen & rh)735 void lp_constraints::get_lr(int index, vecteur &lh, gen &rh) {
736 assert(index<nrows());
737 lh=*lhs[index]._VECTptr;
738 rh=rhs[index];
739 }
740
741 /*
742 * Get the constraint with specified index.
743 */
get(int index,vecteur & lh,gen & rh,int & relation_type)744 void lp_constraints::get(int index, vecteur &lh, gen &rh, int &relation_type) {
745 get_lr(index,lh,rh);
746 relation_type=rv[index];
747 }
748
749 /*
750 * Divide the constraint by g.
751 */
div(int index,const gen & g,GIAC_CONTEXT)752 void lp_constraints::div(int index, const gen &g,GIAC_CONTEXT) {
753 assert(index<nrows() && !is_zero(g));
754 lhs[index]=divvecteur(*lhs[index]._VECTptr,g);
755 rhs[index]=rhs[index]/g;
756 if (is_strictly_positive(-g,contextptr))
757 rv[index]*=-1;
758 }
759
760 /*
761 * Subtract vector from the constraint.
762 */
subtract(int index,const vecteur & v,const gen & g)763 void lp_constraints::subtract(int index, const vecteur &v, const gen &g) {
764 assert(index<nrows());
765 lhs[index]=subvecteur(*lhs[index]._VECTptr,v);
766 rhs[index]-=g;
767 }
768
769 /*
770 * Remove row with specified index.
771 */
remove(int index)772 void lp_constraints::remove(int index) {
773 lhs.erase(lhs.begin()+index);
774 rhs.erase(rhs.begin()+index);
775 rv.erase(rv.begin()+index);
776 if (int(score.size())>index)
777 score.erase(score.begin()+index);
778 }
779
780 /*
781 * Add identifiers from g to variable_identifiers.
782 */
add_identifiers_from(const gen & g)783 void lp_problem::add_identifiers_from(const gen &g) {
784 vecteur vars(*_lname(g,ctx)._VECTptr);
785 for (const_iterateur it=vars.begin();it!=vars.end();++it) {
786 if (!contains(variable_identifiers,*it))
787 variable_identifiers.push_back(*it);
788 }
789 variable_identifiers=*_sort(variable_identifiers,ctx)._VECTptr;
790 }
791
792 /*
793 * Return variable index corresponding to the given identifier.
794 */
get_variable_index(const identificateur & idnt)795 int lp_problem::get_variable_index(const identificateur &idnt) {
796 int n=variable_identifiers.size();
797 for (int i=0;i<n;++i) {
798 if (*variable_identifiers[i]._IDNTptr==idnt)
799 return i;
800 }
801 return -1;
802 }
803
804 /*
805 * Create the problem variables. Assume that they are continuous and
806 * unrestricted by default.
807 */
create_variables(int n)808 void lp_problem::create_variables(int n) {
809 variables=vector<lp_variable>(n);
810 nvars_initial=n;
811 for (int i=0;i<n;++i) {
812 lp_variable var;
813 var.range.lbound=minf; // default: no restrictions whatsoever
814 variables[i]=var;
815 }
816 }
817
818 /*
819 * Tighten both upper and lower bound of the variable.
820 */
tighten_variable_bounds(int i,const gen & l,const gen & u)821 void lp_problem::tighten_variable_bounds(int i, const gen &l, const gen &u) {
822 lp_variable &var=variables[i];
823 var.range.tighten_lbound(l,ctx);
824 var.range.tighten_ubound(u,ctx);
825 }
826
827 /*
828 * Output solution in form [x1=v1,x2=v2,...,xn=vn] where xk are variable
829 * identifiers and vk are solution values.
830 */
output_solution()831 vecteur lp_problem::output_solution() {
832 if (variable_identifiers.empty())
833 return solution;
834 return *_zip(makesequence(at_equal,variable_identifiers,solution),ctx)._VECTptr;
835 }
836
837 /*
838 * Determine coeffcients of linear combination g of variables x in
839 * variable_identifiers. varcoeffs C and freecoeff c are filled such that
840 * g=C*x+c.
841 */
lincomb_coeff(const gen & g,vecteur & varcoeffs,gen & freecoeff)842 bool lp_problem::lincomb_coeff(const gen &g,vecteur &varcoeffs,gen &freecoeff) {
843 gen e(g),a;
844 varcoeffs.clear();
845 for (const_iterateur it=variable_identifiers.begin();it!=variable_identifiers.end();++it) {
846 a=gen(0);
847 if (is_constant_wrt(e,*it,ctx) || (is_linear_wrt(e,*it,a,e,ctx) && is_realcons(a,ctx)))
848 varcoeffs.push_back(a);
849 else return false;
850 }
851 return is_realcons(freecoeff=e,ctx);
852 }
853
854 /*
855 * Add slack variables to the problem if necessary (i.e. convert all
856 * inequalities to equalities).
857 */
add_slack_variables()858 void lp_problem::add_slack_variables() {
859 ints posv;
860 int nv0=constr.ncols();
861 for (int i=0;i<nc();++i) {
862 if (constr.rv[i]==_LP_EQ)
863 continue;
864 append_column(constr.lhs,multvecteur(gen(-constr.rv[i]),singleton(nc(),i)));
865 constr.rv[i]=_LP_EQ;
866 variables.push_back(lp_variable()); //add slack variable
867 posv.push_back(i);
868 }
869 objective.first.resize(nv());
870 //determine types of slack variables
871 vecteur lh;
872 gen rh;
873 for (int k=0;k<int(posv.size());++k) {
874 int i=posv[k],j=nv0+k;
875 lp_variable &var=variables[j];
876 constr.get_lr(i,lh,rh);
877 if (is_exact(lh) && is_exact(rh)) {
878 gen den(_denom(rh,ctx));
879 for (int l=0;l<nv0;++l) {
880 if (is_zero(lh[l]))
881 continue;
882 if (!variables[l].is_integral) {
883 den=undef;
884 break;
885 }
886 if (!lh[l].is_integer())
887 den=_lcm(makesequence(den,_denom(lh[l],ctx)),ctx);
888 }
889 if (!is_undef(den)) {
890 lh=multvecteur(den,lh);
891 rh=den*rh;
892 var.is_integral=true;
893 }
894 }
895 }
896 }
897
898 /*
899 * Make all decision variables bounded below by negating variables unrestricted
900 * below and replacing variables unrestricted from both above and below with
901 * the difference of two nonnegative variables. This process is reversed after
902 * an optimal solution is found.
903 */
make_all_vars_bounded_below()904 void lp_problem::make_all_vars_bounded_below() {
905 for (int i=nv()-1;i>=0;--i) {
906 lp_variable &var=variables[i];
907 if (var.range.is_unrestricted_below()) {
908 if (var.range.is_unrestricted_above()) {
909 var.range.lbound=gen(0);
910 lp_variable negvar(var);
911 var.sign_type=_LP_VARSIGN_POS_PART;
912 negvar.sign_type=_LP_VARSIGN_NEG_PART;
913 variables.insert(variables.begin()+i,negvar);
914 objective.first.insert(objective.first.begin()+i,-objective.first[i]);
915 constr.duplicate_column(i);
916 }
917 else {
918 var.range.lbound=-var.range.ubound;
919 var.range.ubound=inf;
920 var.sign_type=_LP_VARSIGN_NEG;
921 objective.first[i]=-objective.first[i];
922 }
923 constr.negate_column(i);
924 }
925 }
926 }
927
928 /*
929 * Force all problem parameters to be exact.
930 */
make_problem_exact()931 void lp_problem::make_problem_exact() {
932 objective.first=*exact(objective.first,ctx)._VECTptr;
933 objective.second=exact(objective.second,ctx);
934 constr.lhs=*exact(constr.lhs,ctx)._VECTptr;
935 constr.rhs=*exact(constr.rhs,ctx)._VECTptr;
936 for (vector<lp_variable>::iterator it=variables.begin();it!=variables.end();++it) {
937 it->range.lbound=exact(it->range.lbound,ctx);
938 it->range.ubound=exact(it->range.ubound,ctx);
939 }
940 }
941
942 /*
943 * Report status.
944 */
report_status(const char * msg,int count)945 void lp_problem::report_status(const char *msg, int count) {
946 char buf[16];
947 sprintf(buf,"%d: ",count);
948 int nd=numdigits((unsigned)count);
949 string str(msg);
950 str.insert(0,buf);
951 while (nd<8) {
952 str.insert(str.begin(),(char)32);
953 ++nd;
954 }
955 message(str.c_str());
956 }
957
958 /*
959 * Solve the problem using the specified settings.
960 */
solve()961 int lp_problem::solve() {
962 stats=lp_stats();
963 char buffer[1024];
964 make_problem_exact();
965 add_slack_variables();
966 make_all_vars_bounded_below();
967 objective_norm=_evalf(_l2norm(objective.first,ctx),ctx).DOUBLE_val();
968 if (settings.maximize)
969 objective.first=-objective.first;
970 message("Optimizing...");
971 int result;
972 optimum=undef;
973 double opt_approx;
974 lp_node root;
975 root.prob=this;
976 root.ranges=vector<lp_range>(nv());
977 root.depth=0;
978 if ((result=root.solve_relaxation())!=_LP_SOLVED)
979 return result; //optimal solution does not exist
980 if (root.is_integer_feasible()) {
981 solution=root.solution;
982 optimum=root.optimum;
983 }
984 if (is_undef(optimum)) {
985 message("Applying branch&bound method to find integer feasible solutions...");
986 double root_optimum=_evalf(root.optimum,ctx).DOUBLE_val();
987 double root_infeas=_evalf(root.infeas,ctx).DOUBLE_val();
988 root.ranges.resize(nv());
989 vector<lp_node> active_nodes;
990 active_nodes.push_back(root);
991 clock_t t=clock(),t0=t,now;
992 int n,j,k,depth;
993 double opt_lbound,fr,max_score;
994 bool depth_exceeded=false,incumbent_updated,is_use_pseudocost=false;
995 map<double,int> candidates;
996 gen lb,ub;
997 while (!active_nodes.empty()) {
998 if (stats.subproblems_examined>=settings.node_limit) {
999 message("Warning: node limit exceeded",true);
1000 break;
1001 }
1002 n=active_nodes.size();
1003 if (n>stats.max_active_nodes)
1004 stats.max_active_nodes=n;
1005 opt_lbound=DBL_MAX;
1006 k=-1;
1007 depth=-1;
1008 for (int i=0;i<n;++i) {
1009 lp_node &node=active_nodes[i];
1010 if (opt_lbound>node.opt_approx) {
1011 opt_lbound=node.opt_approx;
1012 k=i;
1013 }
1014 if (settings.nodeselect==_LP_BREADTHFIRST ||
1015 (settings.nodeselect==_LP_HYBRID && !is_undef(optimum)) ||
1016 node.depth<depth)
1017 continue;
1018 if (node.depth>depth) {
1019 candidates.clear();
1020 depth=node.depth;
1021 }
1022 candidates[node.opt_approx]=i;
1023 }
1024 if (settings.nodeselect==_LP_DEPTHFIRST ||
1025 (settings.nodeselect==_LP_HYBRID && is_undef(optimum)))
1026 k=candidates.begin()->second;
1027 if (settings.nodeselect==_LP_BEST_PROJECTION) {
1028 double bestproj=DBL_MAX,proj,iopt=is_undef(optimum)?0:_evalf(optimum,ctx).DOUBLE_val();
1029 for (int i=0;i<n;++i) {
1030 lp_node &node=active_nodes[i];
1031 proj=_evalf(node.optimum,ctx).DOUBLE_val()+
1032 (iopt-root_optimum)*_evalf(node.infeas,ctx).DOUBLE_val()/root_infeas;
1033 if (proj<bestproj) {
1034 bestproj=proj;
1035 k=i;
1036 }
1037 }
1038 }
1039 if (k<0) {
1040 message("Error: node selection strategy failed",true);
1041 break;
1042 }
1043 j=-1;
1044 if (settings.varselect==_LP_PSEUDOCOST || settings.varselect<0) {
1045 max_score=0;
1046 for (int i=0;i<nv();++i) {
1047 if (!active_nodes[k].is_var_fractional(i))
1048 continue;
1049 lp_variable &var=variables[i];
1050 double score=var.score(active_nodes[k].fractional_vars[i]);
1051 if (score==0) {
1052 j=-1;
1053 break;
1054 }
1055 if (score>max_score) {
1056 j=i;
1057 max_score=score;
1058 }
1059 }
1060 if (j>=0 && !is_use_pseudocost) {
1061 report_status("Switched to pseudocost based branching",stats.subproblems_examined);
1062 is_use_pseudocost=true;
1063 }
1064 }
1065 if (j<0) {
1066 switch (settings.varselect) {
1067 case -1:
1068 case _LP_MOSTFRACTIONAL:
1069 case _LP_PSEUDOCOST:
1070 j=active_nodes[k].most_fractional;
1071 break;
1072 case _LP_FIRSTFRACTIONAL:
1073 j=active_nodes[k].fractional_vars.begin()->first;
1074 break;
1075 case _LP_LASTFRACTIONAL:
1076 j=active_nodes[k].fractional_vars.rbegin()->first;
1077 break;
1078 }
1079 }
1080 if (j<0) {
1081 message("Error: branching variable selection strategy failed",true);
1082 break;
1083 }
1084 if (!is_undef(optimum)) {
1085 stats.mip_gap=is_zero(optimum)?-opt_lbound:(opt_approx-opt_lbound)/std::abs(opt_approx);
1086 if (stats.mip_gap<=settings.relative_gap_tolerance) {
1087 if (settings.relative_gap_tolerance>0)
1088 message("Warning: integrality gap threshold reached",true);
1089 break;
1090 }
1091 }
1092 fr=_evalf(active_nodes[k].solution[j],ctx).DOUBLE_val();
1093 lb=gen(int(std::ceil(fr)));
1094 ub=gen(int(std::floor(fr)));
1095 incumbent_updated=false;
1096 lp_node child_node;
1097 for (int i=0;i<2;++i) {
1098 child_node=active_nodes[k].create_child();
1099 if (child_node.depth>settings.depth_limit) {
1100 if (!depth_exceeded) {
1101 message ("Warning: depth limit exceeded",true);
1102 depth_exceeded=true;
1103 }
1104 break;
1105 }
1106 lp_range &range=child_node.ranges[j];
1107 if (i==1) {
1108 range.tighten_lbound(lb,ctx);
1109 if (is_strictly_positive(child_node.ranges[j].lbound,ctx)) {
1110 switch (variables[j].sign_type) {
1111 case _LP_VARSIGN_POS_PART:
1112 child_node.ranges[j-1].ubound=gen(0);
1113 break;
1114 case _LP_VARSIGN_NEG_PART:
1115 child_node.ranges[j+1].ubound=gen(0);
1116 }
1117 }
1118 }
1119 else
1120 range.tighten_ubound(ub,ctx);
1121 ++stats.subproblems_examined;
1122 if (child_node.solve_relaxation()==_LP_SOLVED) {
1123 double p=child_node.fractional_vars[j];
1124 variables[j].update_pseudocost(std::abs(child_node.opt_approx-active_nodes[k].opt_approx),
1125 i==0?p:1-p,i);
1126 if (!is_undef(optimum) && is_greater(child_node.optimum,optimum,ctx))
1127 continue;
1128 if (child_node.is_integer_feasible()) {
1129 //new potential incumbent found
1130 if (is_undef(optimum) || is_strictly_greater(optimum,child_node.optimum,ctx)) {
1131 if (is_undef(optimum))
1132 report_status("Incumbent solution found",stats.subproblems_examined);
1133 else {
1134 sprintf(buffer,"Incumbent solution updated, objective value improvement: %g%%",
1135 (opt_approx-child_node.opt_approx)/std::abs(opt_approx)*100.0);
1136 report_status(buffer,stats.subproblems_examined);
1137 }
1138 solution=child_node.solution;
1139 optimum=child_node.optimum;
1140 opt_approx=child_node.opt_approx;
1141 incumbent_updated=true;
1142 }
1143 }
1144 else
1145 active_nodes.push_back(child_node);
1146 }
1147 }
1148 //fathom
1149 active_nodes.erase(active_nodes.begin()+k);
1150 if (incumbent_updated) {
1151 for (int i=int(active_nodes.size())-1;i>=0;--i) {
1152 if (is_greater(active_nodes[i].optimum,optimum,ctx))
1153 active_nodes.erase(active_nodes.begin()+i);
1154 }
1155 }
1156 now=clock();
1157 if (1e3*double(now-t0)/CLOCKS_PER_SEC>settings.time_limit) {
1158 message("Warning: time limit exceeded",true);
1159 break;
1160 }
1161 if (CLOCKS_PER_SEC/double(now-t)<=settings.status_report_freq) { //report status
1162 sprintf(buffer,"%d nodes active, lower bound: %g",(int)active_nodes.size(),opt_lbound);
1163 string str(buffer);
1164 if (stats.mip_gap>=0) {
1165 sprintf(buffer,", integrality gap: %g%%",stats.mip_gap*100);
1166 str+=string(buffer);
1167 }
1168 report_status(str.c_str(),stats.subproblems_examined);
1169 t=clock();
1170 }
1171 }
1172 //show branch&bound summary
1173 sprintf(buffer,"Summary:\n * %d subproblem(s) examined\n * max. tree size: %d nodes\n * %d Gomory cut(s) applied",
1174 stats.subproblems_examined,stats.max_active_nodes,stats.cuts_applied);
1175 message(buffer);
1176 }
1177 if (is_undef(optimum))
1178 return _LP_INFEASIBLE;
1179 if (settings.maximize)
1180 optimum=-optimum;
1181 optimum+=objective.second;
1182 for (int i=nv()-1;i>=0;--i) {
1183 lp_variable &var=variables[i];
1184 switch (var.sign_type) {
1185 case _LP_VARSIGN_NEG:
1186 solution[i]=-solution[i];
1187 break;
1188 case _LP_VARSIGN_NEG_PART:
1189 solution[i+1]-=solution[i];
1190 solution.erase(solution.begin()+i);
1191 break;
1192 }
1193 }
1194 solution.resize(nvars_initial);
1195 return _LP_SOLVED;
1196 }
1197
1198 #ifdef HAVE_LIBGLPK
1199
1200 /*
1201 * Create GLPK problem from constr.
1202 */
glpk_initialize()1203 glp_prob *lp_problem::glpk_initialize() {
1204 glp_prob *glp=glp_create_prob();
1205 glp_add_rows(glp,nc());
1206 glp_add_cols(glp,nv());
1207 glp_set_obj_dir(glp,settings.maximize?GLP_MAX:GLP_MIN);
1208 for (int i=0;i<=nv();++i) {
1209 glp_set_obj_coef(glp,i,gen2double(i==0?objective.second:objective.first[i-1],ctx));
1210 if (i>0) {
1211 lp_variable &var=variables[i-1];
1212 glp_set_col_kind(glp,i,var.is_integral?GLP_IV:GLP_CV);
1213 int bound_type=GLP_FR;
1214 if (!var.range.is_unrestricted_below() && var.range.is_unrestricted_above())
1215 bound_type=GLP_LO;
1216 else if (var.range.is_unrestricted_below() && !var.range.is_unrestricted_above())
1217 bound_type=GLP_UP;
1218 else if (!var.range.is_unrestricted_below() && !var.range.is_unrestricted_above())
1219 bound_type=GLP_DB;
1220 else if (is_zero(var.range.ubound-var.range.lbound))
1221 bound_type=GLP_FX;
1222 double lo=var.range.is_unrestricted_below()?0.0:gen2double(var.range.lbound,ctx);
1223 double hi=var.range.is_unrestricted_above()?0.0:gen2double(var.range.ubound,ctx);
1224 glp_set_col_bnds(glp,i,bound_type,lo,hi);
1225 }
1226 }
1227 int n=constr.nrows()*constr.ncols();
1228 int *ia=new int[n+1],*ja=new int[n+1]; int k=0;
1229 double *ar=new double[n+1];
1230 gen a;
1231 for (int i=0;i<constr.nrows();++i) {
1232 for (int j=0;j<constr.ncols();++j) {
1233 if (is_zero(a=constr.lhs[i][j]))
1234 continue;
1235 ++k;
1236 ia[k]=i+1;
1237 ja[k]=j+1;
1238 ar[k]=gen2double(a,ctx);
1239 }
1240 double rh=gen2double(constr.rhs[i],ctx);
1241 switch (constr.rv[i]) {
1242 case _LP_EQ:
1243 glp_set_row_bnds(glp,i+1,GLP_FX,rh,0.0);
1244 break;
1245 case _LP_LEQ:
1246 glp_set_row_bnds(glp,i+1,GLP_UP,0.0,rh);
1247 break;
1248 case _LP_GEQ:
1249 glp_set_row_bnds(glp,i+1,GLP_LO,rh,0.0);
1250 break;
1251 }
1252 }
1253 glp_load_matrix(glp,k,ia,ja,ar);
1254 delete [] ia;
1255 delete [] ja;
1256 delete [] ar;
1257 return glp;
1258 }
1259
1260 /*
1261 * Solve LP problem using GLPK implementation of simplex method.
1262 */
glpk_simplex(glp_prob * prob)1263 int lp_problem::glpk_simplex(glp_prob *prob) {
1264 glp_smcp parm;
1265 glp_init_smcp(&parm);
1266 parm.msg_lev=settings.verbose?GLP_MSG_ALL:GLP_MSG_ERR;
1267 parm.it_lim=settings.iteration_limit;
1268 parm.tm_lim=settings.time_limit;
1269 parm.presolve=settings.presolve?GLP_ON:GLP_OFF;
1270 return glp_simplex(prob,&parm);
1271 }
1272
1273 /*
1274 * Solve LP problem using GLPK implementation of interior point method.
1275 */
glpk_interior_point(glp_prob * prob)1276 int lp_problem::glpk_interior_point(glp_prob *prob) {
1277 glp_iptcp parm;
1278 glp_init_iptcp(&parm);
1279 parm.msg_lev=settings.verbose?GLP_MSG_ALL:GLP_MSG_ERR;
1280 return glp_interior(prob,&parm);
1281 }
1282
1283 /*
1284 * Solve MIP problem using GLPK implementation of branch and cut method.
1285 */
glpk_branchcut(glp_prob * prob)1286 int lp_problem::glpk_branchcut(glp_prob *prob) {
1287 glp_iocp parm;
1288 glp_init_iocp(&parm);
1289 parm.msg_lev=settings.verbose?GLP_MSG_ALL:GLP_MSG_ERR;
1290 parm.tm_lim=settings.time_limit;
1291 parm.out_frq=(int)(1000.0/settings.status_report_freq);
1292 parm.mip_gap=settings.relative_gap_tolerance;
1293 parm.gmi_cuts=settings.max_cuts>0?GLP_ON:GLP_OFF;
1294 parm.mir_cuts=settings.max_cuts>0?GLP_ON:GLP_OFF;
1295 parm.clq_cuts=settings.has_binary_vars?GLP_ON:GLP_OFF;
1296 parm.cov_cuts=settings.has_binary_vars?GLP_ON:GLP_OFF;
1297 parm.presolve=GLP_ON;
1298 switch (settings.varselect) {
1299 case _LP_FIRSTFRACTIONAL:
1300 parm.br_tech=GLP_BR_FFV;
1301 break;
1302 case _LP_LASTFRACTIONAL:
1303 parm.br_tech=GLP_BR_LFV;
1304 break;
1305 case _LP_MOSTFRACTIONAL:
1306 parm.br_tech=GLP_BR_MFV;
1307 break;
1308 case _LP_PSEUDOCOST:
1309 parm.br_tech=GLP_BR_PCH;
1310 break;
1311 default:
1312 parm.br_tech=GLP_BR_DTH;
1313 }
1314 switch (settings.nodeselect) {
1315 case _LP_DEPTHFIRST:
1316 parm.bt_tech=GLP_BT_DFS;
1317 break;
1318 case _LP_BREADTHFIRST:
1319 parm.bt_tech=GLP_BT_BFS;
1320 break;
1321 case _LP_BEST_PROJECTION:
1322 parm.bt_tech=GLP_BT_BPH;
1323 break;
1324 default:
1325 parm.bt_tech=GLP_BT_BLB;
1326 }
1327 return glp_intopt(prob,&parm);
1328 }
1329
1330 #endif
1331
1332 /*
1333 * Solve the problem using the GLPK library.
1334 */
glpk_solve()1335 int lp_problem::glpk_solve() {
1336 #ifndef HAVE_LIBGLPK
1337 message("Warning: solving in floating-point arithmetic requires GLPK library",true);
1338 return solve();
1339 #endif
1340 #ifdef HAVE_LIBGLPK
1341 glp_prob *prob=glpk_initialize();
1342 int result,solution_status;
1343 bool is_mip=has_integral_variables();
1344 switch (settings.solver) {
1345 case _LP_SIMPLEX:
1346 if (is_mip) {
1347 result=glpk_branchcut(prob);
1348 if (result==GLP_EMIPGAP)
1349 result=0;
1350 if (result==0)
1351 solution_status=glp_mip_status(prob);
1352 }
1353 else if (!is_mip && (result=glpk_simplex(prob))==0)
1354 solution_status=glp_get_status(prob);
1355 break;
1356 case _LP_INTERIOR_POINT:
1357 if ((result=glpk_interior_point(prob))==0)
1358 solution_status=glp_ipt_status(prob);
1359 break;
1360 }
1361 if (result==0) { //solving process finished successfully
1362 switch (solution_status) {
1363 case GLP_OPT:
1364 solution_status=_LP_SOLVED;
1365 break;
1366 case GLP_FEAS:
1367 message("Warning: solution is possibly not optimal",true);
1368 solution_status=_LP_SOLVED;
1369 break;
1370 case GLP_INFEAS:
1371 message("Warning: solution is infeasible",true);
1372 solution_status=_LP_SOLVED;
1373 break;
1374 case GLP_NOFEAS:
1375 solution_status=_LP_INFEASIBLE;
1376 break;
1377 case GLP_UNBND:
1378 solution_status=_LP_UNBOUNDED;
1379 break;
1380 case GLP_UNDEF:
1381 solution_status=_LP_ERROR;
1382 break;
1383 }
1384 if (solution_status==_LP_SOLVED) { //get solution and optimum
1385 solution=vecteur(nv());
1386 switch (settings.solver) {
1387 case _LP_SIMPLEX:
1388 if (is_mip) {
1389 optimum=gen(glp_mip_obj_val(prob));
1390 for (int i=0;i<nv();++i) {
1391 solution[i]=gen(glp_mip_col_val(prob,i+1));
1392 }
1393 }
1394 else {
1395 optimum=gen(glp_get_obj_val(prob));
1396 for (int i=0;i<nv();++i) {
1397 solution[i]=gen(glp_get_col_prim(prob,i+1));
1398 }
1399 }
1400 break;
1401 case _LP_INTERIOR_POINT:
1402 optimum=gen(glp_ipt_obj_val(prob));
1403 for (int i=0;i<nv();++i) {
1404 solution[i]=gen(glp_ipt_col_prim(prob,i+1));
1405 }
1406 break;
1407 }
1408 }
1409 }
1410 glp_delete_prob(prob);
1411 if (result==0)
1412 return solution_status;
1413 return _LP_ERROR;
1414 #endif
1415 }
1416
1417 /*
1418 * Load problem from file in LP or (possibly gzipped) MPS format. File types
1419 * are distinguished by examining the filename extension. If it is .gz or .mps,
1420 * it is assumed that the file is in (gzipped) MPS format. Else, it is assumed
1421 * that file is in LP format.
1422 */
glpk_load_from_file(const char * fname)1423 bool lp_problem::glpk_load_from_file(const char *fname) {
1424 #ifndef HAVE_LIBGLPK
1425 message("Error: loading CPLEX and MPS files requires GLPK library",true);
1426 return false;
1427 #endif
1428 #ifdef HAVE_LIBGLPK
1429 glp_prob *prob=glp_create_prob();
1430 int result;
1431 string str(fname);
1432 string ext(str.substr(str.find_last_of(".") + 1));
1433 if (ext=="mps" || ext=="gz") //MPS format (modern)
1434 result=glp_read_mps(prob,GLP_MPS_FILE,NULL,fname);
1435 else //LP format
1436 result=glp_read_lp(prob,NULL,fname);
1437 if (result==0) { //successfully loaded
1438 char buffer[1024];
1439 int obj_dir=glp_get_obj_dir(prob),len,t,k;
1440 int nr=glp_get_num_rows(prob),n=glp_get_num_cols(prob);
1441 if (nr*n>LP_CONSTR_MAXSIZE) {
1442 message("Error: constraint matrix is too big to be loaded. Aborting...",true);
1443 return false;
1444 }
1445 int *ind=new int[1+n];
1446 double *val=new double[1+n];
1447 create_variables(n);
1448 constr.lhs=*_matrix(makesequence(gen(nr),gen(nv())),ctx)._VECTptr;
1449 constr.rhs=vecteur(nr);
1450 constr.rv=ints(nr);
1451 variable_identifiers=vecteur(nv());
1452 settings.maximize=(obj_dir==GLP_MAX);
1453 objective.second=double2fraction(glp_get_obj_coef(prob,0),ctx);
1454 objective.first=vecteur(nv());
1455 obj_approx=vector<double>(nv());
1456 for (int j=0;j<nv();++j) {
1457 objective.first[j]=double2fraction(glp_get_obj_coef(prob,j+1),ctx);
1458 obj_approx[j]=std::abs(glp_get_obj_coef(prob,j+1));
1459 lp_variable &var=variables[j];
1460 t=glp_get_col_type(prob,j+1);
1461 k=glp_get_col_kind(prob,j+1);
1462 var.name=string(glp_get_col_name(prob,j+1));
1463 variable_identifiers[j]=identificateur(var.name);
1464 if (t!=GLP_UP && t!=GLP_FR)
1465 var.range.lbound=double2fraction(glp_get_col_lb(prob,j+1),ctx);
1466 if (t!=GLP_LO && t!=GLP_FR)
1467 var.range.ubound=double2fraction(glp_get_col_ub(prob,j+1),ctx);
1468 if (k!=GLP_CV)
1469 var.is_integral=true;
1470 if (k==GLP_BV) {
1471 var.range.tighten_lbound(gen(0),ctx);
1472 var.range.tighten_ubound(gen(1),ctx);
1473 }
1474 }
1475 for (int i=nr;i>0;--i) {
1476 if ((len=glp_get_mat_row(prob,i,ind,val))>0) {
1477 vecteur &row=*constr.lhs[i-1]._VECTptr;
1478 for (int j=1;j<=len;++j) {
1479 row[ind[j]-1]=double2fraction(val[j],ctx);
1480 }
1481 }
1482 t=glp_get_row_type(prob,i);
1483 if (t==GLP_FR) {
1484 sprintf(buffer,"Warning: auxiliary variable bounds not set, discarding constraint %d",i);
1485 message(buffer,true);
1486 constr.remove(i-1);
1487 continue;
1488 }
1489 if (t==GLP_FX) {
1490 constr.rhs[i-1]=double2fraction(glp_get_row_lb(prob,i),ctx);
1491 constr.rv[i-1]=_LP_EQ;
1492 }
1493 else {
1494 if (t==GLP_LO || t==GLP_DB) {
1495 constr.rhs[i-1]=double2fraction(glp_get_row_lb(prob,i),ctx);
1496 constr.rv[i-1]=_LP_GEQ;
1497 }
1498 if (t==GLP_DB)
1499 constr.append(*constr.lhs[i-1]._VECTptr,double2fraction(glp_get_row_ub(prob,i),ctx),_LP_LEQ);
1500 if (t==GLP_UP || t==GLP_DB) {
1501 constr.rhs[i-1]=double2fraction(glp_get_row_ub(prob,i),ctx);
1502 constr.rv[i-1]=_LP_LEQ;
1503 }
1504 }
1505 }
1506 delete [] ind;
1507 delete [] val;
1508 stringstream ss;
1509 const char* name=glp_get_prob_name(prob);
1510 if (name==NULL)
1511 ss << "a problem";
1512 else
1513 ss << "problem \"" << name << "\"";
1514 sprintf(buffer,"Successfully loaded %s with %d variables and %d constraints",
1515 ss.str().c_str(),nv(),nc());
1516 message(buffer);
1517 }
1518 glp_delete_prob(prob);
1519 return !(bool)result;
1520 #endif
1521 }
1522
set_type(int t,GIAC_CONTEXT)1523 void lp_variable::set_type(int t,GIAC_CONTEXT) {
1524 switch (t) {
1525 case _LP_BINARYVARIABLES:
1526 range.tighten_lbound(gen(0),contextptr);
1527 range.tighten_ubound(gen(1),contextptr);
1528 case _LP_INTEGERVARIABLES:
1529 is_integral=true;
1530 break;
1531 }
1532 }
1533
assign_variable_types(const gen & g,int t,lp_problem & prob)1534 bool assign_variable_types(const gen &g,int t,lp_problem &prob) {
1535 pair<gen,gen> range;
1536 int i,i0=array_start(prob.ctx);
1537 if (g.type==_VECT) {
1538 for (const_iterateur it=g._VECTptr->begin();it!=g._VECTptr->end();++it) {
1539 if (!assign_variable_types(*it,t,prob))
1540 return false;
1541 }
1542 }
1543 else if ((g.type==_IDNT && (i=prob.get_variable_index(*g._IDNTptr))>=0) ||
1544 (g.is_integer() && (i=g.val-i0)>=0))
1545 prob.variables[i].set_type(t,prob.ctx);
1546 else if (interval2pair(g,range,prob.ctx) &&
1547 range.first.is_integer() && range.second.is_integer()) {
1548 for (i=range.first.val;i<=range.second.val;++i) {
1549 prob.variables[i-i0].set_type(t,prob.ctx);
1550 }
1551 }
1552 else
1553 return false;
1554 return true;
1555 }
1556
parse_limit(const gen & g,int & lim,GIAC_CONTEXT)1557 void parse_limit(const gen &g,int &lim,GIAC_CONTEXT) {
1558 if (is_positive(g,contextptr)) {
1559 if (g.is_integer())
1560 lim=g.val;
1561 else if (is_inf(g))
1562 lim=RAND_MAX;
1563 }
1564 }
1565
parse_options_and_bounds(const_iterateur & it,const_iterateur & itend,lp_problem & prob)1566 bool parse_options_and_bounds(const_iterateur &it,const_iterateur &itend,lp_problem &prob) {
1567 for (;it!=itend;++it) {
1568 if (*it==at_maximize)
1569 prob.settings.maximize=true;
1570 else if (it->is_integer()) {
1571 switch(it->val) {
1572 case _LP_MAXIMIZE:
1573 prob.settings.maximize=true;
1574 break;
1575 case _LP_VERBOSE:
1576 prob.settings.verbose=true;
1577 break;
1578 }
1579 }
1580 else if (it->is_symb_of_sommet(at_equal)) {
1581 //parse the argument in form "option=value"
1582 pair<gen,gen> bounds;
1583 gen a,lh(_lhs(*it,prob.ctx)),rh(_rhs(*it,prob.ctx));
1584 if (lh.type==_IDNT && interval2pair(rh,bounds,prob.ctx)) {
1585 int vi=prob.get_variable_index(*lh._IDNTptr);
1586 if (vi>=0)
1587 prob.tighten_variable_bounds(vi,bounds.first,bounds.second);
1588 }
1589 else if (lh==at_maximize && rh.is_integer())
1590 prob.settings.maximize=(bool)rh.to_int();
1591 else if (lh==at_assume && rh.is_integer())
1592 prob.settings.assumption=rh.val;
1593 else if (lh.is_integer()) {
1594 switch(lh.val) {
1595 case _LP_INTEGERVARIABLES:
1596 case _LP_BINARYVARIABLES:
1597 if (!assign_variable_types(rh,lh.val,prob))
1598 return false;
1599 break;
1600 case _LP_DEPTHLIMIT:
1601 parse_limit(rh,prob.settings.depth_limit,prob.ctx);
1602 break;
1603 case _LP_NODE_LIMIT:
1604 parse_limit(rh,prob.settings.node_limit,prob.ctx);
1605 break;
1606 case _LP_ITERATION_LIMIT:
1607 parse_limit(rh,prob.settings.iteration_limit,prob.ctx);
1608 break;
1609 case _LP_MAX_CUTS:
1610 parse_limit(rh,prob.settings.max_cuts,prob.ctx);
1611 break;
1612 case _LP_TIME_LIMIT:
1613 if (is_strictly_positive(rh,prob.ctx)) {
1614 if (is_inf(rh))
1615 prob.settings.time_limit=RAND_MAX;
1616 a=_evalf(rh,prob.ctx);
1617 if (a.type==_DOUBLE_)
1618 prob.settings.time_limit=int(a.DOUBLE_val()*1000.0);
1619 }
1620 break;
1621 case _LP_NODESELECT:
1622 if (rh.is_integer())
1623 prob.settings.nodeselect=rh.val;
1624 break;
1625 case _LP_VARSELECT:
1626 if (rh.is_integer())
1627 prob.settings.varselect=rh.val;
1628 break;
1629 case _LP_METHOD:
1630 if (rh==at_exact)
1631 prob.settings.precision=_LP_EXACT;
1632 else if (rh==at_float)
1633 prob.settings.precision=_LP_INEXACT;
1634 else if (rh.is_integer() && rh.val==_LP_INTERIOR_POINT) {
1635 prob.settings.precision=_LP_INEXACT;
1636 prob.settings.solver=_LP_INTERIOR_POINT;
1637 }
1638 break;
1639 case _LP_GAP_TOLERANCE:
1640 a=_evalf(rh,prob.ctx);
1641 if (a.type==_DOUBLE_ && is_strictly_positive(a,prob.ctx))
1642 prob.settings.relative_gap_tolerance=a.DOUBLE_val();
1643 break;
1644 case _LP_MAXIMIZE:
1645 if (rh.is_integer())
1646 prob.settings.maximize=(bool)rh.to_int();
1647 break;
1648 case _LP_VERBOSE:
1649 if (rh.is_integer())
1650 prob.settings.verbose=(bool)rh.to_int();
1651 break;
1652 }
1653 }
1654 }
1655 }
1656 return true;
1657 }
1658
parse_constraints(bool is_matrix_form,const_iterateur & it,const_iterateur & itend,lp_problem & prob)1659 int parse_constraints(bool is_matrix_form,const_iterateur &it,const_iterateur &itend,lp_problem &prob) {
1660 vecteur cv;
1661 if (is_matrix_form) {
1662 if (it->type!=_VECT)
1663 return _LP_ERR_TYPE;
1664 cv=*it->_VECTptr;
1665 int n=cv.size();
1666 if ((n % 2)!=0)
1667 return _LP_ERR_SIZE;
1668 if (n>0) {
1669 //constraints are given in form [A,b,Aeq,beq] such that
1670 //A.x<=b and/or Aeq.x=beq (Aeq and beq are optional)
1671 if (cv[0].type!=_VECT || cv[1].type!=_VECT)
1672 return _LP_ERR_TYPE;
1673 vecteur &A=*cv[0]._VECTptr,&b=*cv[1]._VECTptr;
1674 int len;
1675 if ((len=A.size())!=int(b.size()))
1676 return _LP_ERR_DIM;
1677 for (int i=0;i<len;++i) {
1678 prob.constr.append(*A[i]._VECTptr,b[i],_LP_LEQ);
1679 }
1680 if (n>2) { //there are equality constraints
1681 if (cv[2].type!=_VECT || cv[3].type!=_VECT)
1682 return _LP_ERR_TYPE;
1683 vecteur &Aeq=*cv[2]._VECTptr,&beq=*cv[3]._VECTptr;
1684 if ((len=Aeq.size())!=int(beq.size()))
1685 return _LP_ERR_DIM;
1686 for (int i=0;i<len;++i) {
1687 prob.constr.append(*Aeq[i]._VECTptr,beq[i],_LP_EQ);
1688 }
1689 }
1690 if (!ckmatrix(prob.constr.lhs))
1691 return _LP_ERR_DIM;
1692 if (prob.nv()==0)
1693 prob.create_variables(prob.constr.ncols());
1694 else if (prob.nv()!=prob.constr.ncols())
1695 return _LP_ERR_DIM;
1696 }
1697 //if the third arg is [bl,bu] then assume bl<=x<=bu
1698 ++it;
1699 if (it!=itend && it->type==_VECT) {
1700 matrice &bounds=*it->_VECTptr;
1701 if (!ckmatrix(bounds) || int(bounds.size())!=2 ||
1702 int(bounds.front()._VECTptr->size())<prob.nv())
1703 return _LP_ERR_DIM;
1704 if (prob.nv()==0)
1705 prob.create_variables(bounds.front()._VECTptr->size());
1706 for (int i=0;i<prob.nv();++i) {
1707 prob.tighten_variable_bounds(i,bounds[0][i],bounds[1][i]);
1708 }
1709 ++it;
1710 }
1711 if (prob.nv()==0)
1712 return _LP_ERR_SIZE;
1713 if (prob.objective.first.empty())
1714 prob.set_objective(vecteur(prob.nv(),gen(0)),gen(0));
1715 }
1716 else { //the problem is given in symbolic form
1717 if (it->type==_VECT) {
1718 cv=*it->_VECTptr;
1719 prob.add_identifiers_from(cv);
1720 ++it;
1721 }
1722 prob.create_variables(prob.variable_identifiers.size());
1723 vecteur sides,c;
1724 pair<gen,gen> bounds;
1725 gen r;
1726 int rel;
1727 for (const_iterateur jt=cv.begin();jt!=cv.end();++jt) {
1728 if (jt->is_symb_of_sommet(at_equal))
1729 rel=_LP_EQ;
1730 else if (jt->is_symb_of_sommet(at_inferieur_egal))
1731 rel=_LP_LEQ;
1732 else if (jt->is_symb_of_sommet(at_superieur_egal))
1733 rel=_LP_GEQ;
1734 else
1735 return _LP_ERR_TYPE;
1736 sides=*jt->_SYMBptr->feuille._VECTptr;
1737 if (jt->is_symb_of_sommet(at_equal) && interval2pair(sides.back(),bounds,prob.ctx)) {
1738 if (!prob.lincomb_coeff(sides.front(),c,r))
1739 return _LP_ERR_TYPE;
1740 prob.constr.append(c,bounds.first-r,_LP_GEQ);
1741 prob.constr.append(c,bounds.second-r,_LP_LEQ);
1742 continue;
1743 }
1744 if (!prob.lincomb_coeff(sides.front()-sides.back(),c,r))
1745 return _LP_ERR_TYPE;
1746 prob.constr.append(c,-r,rel);
1747 }
1748 }
1749 if (!parse_options_and_bounds(it,itend,prob))
1750 return _LP_ERR_TYPE;
1751 return 0;
1752 }
1753
_lpsolve(const gen & args,GIAC_CONTEXT)1754 gen _lpsolve(const gen &args,GIAC_CONTEXT) {
1755 gen g(args);
1756 if (g.type==_STRNG)
1757 g=makesequence(g);
1758 if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
1759 return gensizeerr(contextptr);
1760 vecteur &gv=*g._VECTptr,obj;
1761 const_iterateur it=gv.begin(),itend=gv.end();
1762 lp_problem prob(contextptr); //create LP problem with default settings
1763 bool is_matrix_form=false;
1764 if (it->type==_STRNG) { //problem is given in file
1765 int len=_size(*it,contextptr).val;
1766 string fname(it->_STRNGptr->begin(),it->_STRNGptr->begin()+len);
1767 if (!prob.glpk_load_from_file(fname.c_str()))
1768 return undef;
1769 ++it;
1770 if (it->type==_VECT)
1771 return gentypeerr(contextptr);
1772 }
1773 else { //problem is entered manually
1774 is_matrix_form=it->type==_VECT;
1775 if (is_matrix_form) {
1776 obj=*it->_VECTptr;
1777 if (!obj.empty()) {
1778 prob.set_objective(obj,gen(0));
1779 prob.create_variables(obj.size());
1780 }
1781 }
1782 else {
1783 prob.add_identifiers_from(*it);
1784 gen obj_ct;
1785 if (!prob.lincomb_coeff(*it,obj,obj_ct))
1786 return _LP_ERR_TYPE;
1787 prob.set_objective(obj,obj_ct);
1788 }
1789 ++it;
1790 }
1791 switch (parse_constraints(is_matrix_form,it,itend,prob)) {
1792 case _LP_ERR_SIZE:
1793 return gensizeerr(contextptr);
1794 case _LP_ERR_TYPE:
1795 return gentypeerr(contextptr);
1796 case _LP_ERR_DIM:
1797 return gendimerr(contextptr);
1798 default:
1799 if (prob.nc()==0) {
1800 for (int i=0;i<prob.nv();++i) {
1801 lp_variable &var=prob.variables[i];
1802 vecteur row(prob.nv(),gen(0));
1803 if (!is_inf(var.range.lbound)) {
1804 row[i]=gen(1);
1805 prob.constr.append(row,var.range.lbound,_LP_GEQ);
1806 break;
1807 }
1808 else if (!is_inf(var.range.ubound)) {
1809 row[i]=gen(1);
1810 prob.constr.append(row,var.range.ubound,_LP_LEQ);
1811 break;
1812 }
1813 }
1814 if (prob.nc()==0) {
1815 prob.message("Error: no constraints detected",true);
1816 return gensizeerr(contextptr);
1817 }
1818 }
1819 }
1820 for (int i=0;i<prob.nv();++i) {
1821 lp_variable &var=prob.variables[i];
1822 switch (prob.settings.assumption) {
1823 case _ZINT:
1824 case _LP_INTEGER:
1825 var.set_type(_LP_INTEGERVARIABLES,contextptr);
1826 break;
1827 case _LP_BINARY:
1828 var.set_type(_LP_BINARYVARIABLES,contextptr);
1829 break;
1830 case _NONNEGINT:
1831 case _LP_NONNEGINT:
1832 var.set_type(_LP_INTEGERVARIABLES,contextptr);
1833 case _LP_NONNEGATIVE:
1834 var.range.tighten_lbound(gen(0),contextptr);
1835 break;
1836 }
1837 }
1838 bool is_solver_exact;
1839 switch (prob.settings.precision) {
1840 case _LP_EXACT:
1841 is_solver_exact=true;
1842 break;
1843 case _LP_INEXACT:
1844 is_solver_exact=false;
1845 break;
1846 case _LP_PROB_DEPENDENT:
1847 is_solver_exact=!prob.has_approx_coefficients();
1848 break;
1849 }
1850 vector<lp_variable>::const_iterator vt=prob.variables.begin();
1851 for (;vt!=prob.variables.end();++vt) {
1852 if (vt->is_integral && is_zero(vt->range.lbound) && is_one(vt->range.ubound)) {
1853 prob.settings.has_binary_vars=true;
1854 break;
1855 }
1856 }
1857 //solve the problem
1858 switch (is_solver_exact?prob.solve():prob.glpk_solve()) {
1859 case _LP_INFEASIBLE:
1860 prob.message("Error: problem has no feasible solutions",true);
1861 return vecteur(0);
1862 case _LP_UNBOUNDED:
1863 prob.message("Error: problem has unbounded solution",true);
1864 return makevecteur(prob.settings.maximize?inf:minf,vecteur(0));
1865 case _LP_ERROR:
1866 return undef;
1867 default: //_LP_SOLVED
1868 return gen(makevecteur(_eval(prob.optimum,contextptr),_eval(prob.output_solution(),contextptr)),_LIST__VECT);
1869 }
1870 }
1871 static const char _lpsolve_s[]="lpsolve";
1872 static define_unary_function_eval (__lpsolve,&_lpsolve,_lpsolve_s);
1873 define_unary_function_ptr5(at_lpsolve,alias_at_lpsolve,&__lpsolve,0,true)
1874
1875 #ifndef NO_NAMESPACE_GIAC
1876 }
1877 #endif //ndef NO_NAMESPACE_GIAC
1878