1 #ifndef __LPSOLVE_H
2 #define __LPSOLVE_H
3 
4 #include "config.h"
5 #include "gen.h"
6 #include "unary.h"
7 #ifdef HAVE_LIBGLPK
8 #include <glpk.h>
9 #endif
10 
11 #ifndef NO_NAMESPACE_GIAC
12 namespace giac {
13 #endif //ndef NO_NAMESPACE_GIAC
14 
15 #define LP_SCORE_FACTOR 0.1667
16 #define LP_MIN_AWAY 0.08
17 #define LP_MIN_PARALLELISM 0.86
18 #define LP_MAX_MAGNITUDE 1e6
19 #define LP_CONSTR_MAXSIZE 1e5
20 
21 typedef std::vector<int> ints;
22 
23 enum lp_results {
24     _LP_SOLVED,
25     _LP_INFEASIBLE,
26     _LP_UNBOUNDED,
27     _LP_ERROR
28 };
29 
30 enum lp_parsing_errors {
31     _LP_ERR_SIZE = 1,
32     _LP_ERR_TYPE = 2,
33     _LP_ERR_DIM  = 3
34 };
35 
36 enum lp_precision_types {
37     _LP_EXACT,
38     _LP_INEXACT,
39     _LP_PROB_DEPENDENT
40 };
41 
42 enum lp_relation_types {
43     _LP_LEQ = -1,
44     _LP_EQ  =  0,
45     _LP_GEQ =  1
46 };
47 
48 enum lp_variable_sign_types {
49     _LP_VARSIGN_POS,
50     _LP_VARSIGN_NEG,
51     _LP_VARSIGN_POS_PART,
52     _LP_VARSIGN_NEG_PART
53 };
54 
55 struct lp_settings {
56     //solver parameters
57     int solver;
58     int precision;
59     bool presolve;
60     bool maximize;
61     int assumption;
62     int iteration_limit;
63     bool has_binary_vars;
64     //branch&bound parameters
65     int varselect;
66     int nodeselect;
67     double relative_gap_tolerance;
68     int depth_limit;
69     int node_limit;
70     int time_limit; //in miliseconds
71     int max_cuts;
72     //message report parameters
73     bool verbose;
74     double status_report_freq;
75     lp_settings();
76 };
77 
78 struct lp_stats {
79     int subproblems_examined;
80     int cuts_applied;
81     int max_active_nodes;
82     double mip_gap;
83     lp_stats();
84 };
85 
86 struct lp_range {
87     gen lbound;
88     gen ubound;
tighten_lboundlp_range89     inline void tighten_lbound(const gen &l,GIAC_CONTEXT) { lbound=max(lbound,l,contextptr); }
tighten_uboundlp_range90     inline void tighten_ubound(const gen &u,GIAC_CONTEXT) { ubound=min(ubound,u,contextptr); }
is_unrestricted_belowlp_range91     inline bool is_unrestricted_below() { return is_inf(lbound); }
is_unrestricted_abovelp_range92     inline bool is_unrestricted_above() { return is_inf(ubound); }
93     lp_range();
94 };
95 
96 struct lp_variable {
97     bool is_integral;
98     int sign_type;
99     lp_range range;
100     std::string name;
101     double pseudocost[2];
102     int nbranch[2];
103     lp_variable();
104     void set_type(int t,GIAC_CONTEXT);
105     void update_pseudocost(double delta,double fr,int dir);
106     double score(double fr);
107 };
108 
109 struct lp_constraints {
110     matrice lhs;
111     vecteur rhs;
112     ints rv;
113     std::vector<double> score;
nrowslp_constraints114     inline int nrows() { return lhs.size(); }
ncolslp_constraints115     inline int ncols() { return lhs.empty()?0:lhs.front()._VECTptr->size(); }
116     void append(const vecteur &lh,const gen &rh,int relation_type);
117     vecteur column(int index);
118     void duplicate_column(int index);
119     void negate_column(int index);
120     void subtract_from_rhs_column(const vecteur &v);
121     void set(int index,const vecteur &lh,const gen &rh,int relation_type);
122     void get(int index,vecteur &lh,gen &rh,int &relation_type);
123     void get_lr(int index,vecteur &lh,gen &rh);
124     void div(int index,const gen &g,GIAC_CONTEXT);
125     void subtract(int index,const vecteur &v,const gen &g);
126     void remove(int index);
127 };
128 
129 struct lp_node;
130 
131 struct lp_problem {
132     const context *ctx;
133     std::pair<vecteur,gen> objective;
134     double objective_norm;
135     std::vector<double> obj_approx;
136     std::vector<lp_variable> variables;
137     vecteur variable_identifiers;
138     lp_constraints constr;
139     lp_constraints cuts;
140     lp_settings settings;
141     lp_stats stats;
142     int nvars_initial;
143     vecteur solution;
144     gen optimum;
lp_problemlp_problem145     lp_problem(GIAC_CONTEXT) {
146         ctx=contextptr;
147         settings=lp_settings();
148     }
nclp_problem149     inline int nc() { return constr.lhs.size(); }
nvlp_problem150     inline int nv() { return variables.size(); }
151     void message(const char* msg,bool err=false);
152     void report_status(const char* msg,int count);
153     void add_identifiers_from(const gen &g);
154     int get_variable_index(const identificateur &idnt);
155     void set_objective(const vecteur &v,const gen &ft);
156     void create_variables(int n);
157     void make_problem_exact();
158     void add_slack_variables();
159     void tighten_variable_bounds(int i,const gen &l,const gen &u);
160     void make_all_vars_bounded_below();
161     bool has_integral_variables();
162     bool has_approx_coefficients();
163     bool lincomb_coeff(const gen &g,vecteur &varcoeffs,gen &freecoeff);
164     int solve();
165     vecteur output_solution();
166     //GLPK routines
167 #ifdef HAVE_LIBGLPK
168     glp_prob *glpk_initialize();
169     int glpk_simplex(glp_prob *prob);
170     int glpk_interior_point(glp_prob *prob);
171     int glpk_branchcut(glp_prob *prob);
172 #endif
173     int glpk_solve();
174     bool glpk_load_from_file(const char *fname);
175 };
176 
177 struct lp_node {
178     lp_problem *prob;
179     int depth;
180     std::vector<lp_range> ranges;
181     gen optimum;
182     vecteur solution;
183     double opt_approx;
184     gen infeas;
185     int most_fractional;
186     std::map<int,double> fractional_vars;
187     ints cut_indices;
is_integer_feasiblelp_node188     inline bool is_integer_feasible() { return is_zero(infeas); }
189     bool is_var_fractional(int index);
190     gen fracpart(const gen &g);
191     lp_node create_child();
192     int solve_relaxation();
193 };
194 
195 gen _lpsolve(const gen &args,GIAC_CONTEXT);
196 extern const unary_function_ptr * const  at_lpsolve;
197 
198 #ifndef NO_NAMESPACE_GIAC
199 } //namespace giac
200 #endif //ndef NO_NAMESPACE_GIAC
201 #endif //__LPSOLVE_H
202