1 /* NUMlinprog.cpp
2 *
3 * Copyright (C) 2008,2011,2012,2015-2020 Paul Boersma
4 *
5 * This code is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This code is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 * See the GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this work. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include "melder.h"
20 #include "../external/glpk/glpk.h"
21
22 struct structNUMlinprog {
23 glp_prob *linearProgram;
24 integer numberOfConstraints, ivar, numberOfVariables;
25 autovector<int> ind;
26 autoVEC val;
27 int status;
28 };
29
NUMlinprog_delete(NUMlinprog me)30 void NUMlinprog_delete (NUMlinprog me) {
31 if (! me)
32 return;
33 if (my linearProgram)
34 glp_delete_prob (my linearProgram);
35 my ind.reset();
36 my val.reset();
37 Melder_free (me);
38 }
39
NUMlinprog_new(bool maximize)40 NUMlinprog NUMlinprog_new (bool maximize) {
41 NUMlinprog me = nullptr;
42 try {
43 me = Melder_calloc (structNUMlinprog, 1);
44 my linearProgram = glp_create_prob (); // TODO: check
45 glp_set_obj_dir (my linearProgram, maximize ? GLP_MAX : GLP_MIN);
46 } catch (MelderError) {
47 if (me) NUMlinprog_delete (me);
48 return NULL;
49 }
50 return me;
51 }
52
NUMlinprog_addVariable(NUMlinprog me,double lowerBound,double upperBound,double coeff)53 void NUMlinprog_addVariable (NUMlinprog me, double lowerBound, double upperBound, double coeff) {
54 glp_add_cols (my linearProgram, 1);
55 glp_set_col_bnds (my linearProgram, (int) ++ my numberOfVariables,
56 isundef (lowerBound) ? ( isundef (upperBound) ? GLP_FR : GLP_UP ) :
57 isundef (upperBound) ? GLP_LO :
58 lowerBound == upperBound ? GLP_FX : GLP_DB,
59 lowerBound, upperBound);
60 glp_set_obj_coef (my linearProgram, (int) my ivar, coeff);
61 }
62
NUMlinprog_addConstraint(NUMlinprog me,double lowerBound,double upperBound)63 void NUMlinprog_addConstraint (NUMlinprog me, double lowerBound, double upperBound) {
64 try {
65 if (NUMisEmpty (my ind.get())) {
66 my ind = newvectorzero<int> (my numberOfVariables);
67 my val = zero_VEC (my numberOfVariables);
68 }
69 glp_add_rows (my linearProgram, 1); // TODO: check
70 glp_set_row_bnds (my linearProgram, (int) ++ my numberOfConstraints,
71 isundef (lowerBound) ? ( isundef (upperBound) ? GLP_FR : GLP_UP ) :
72 isundef (upperBound) ? GLP_LO :
73 lowerBound == upperBound ? GLP_FX : GLP_DB, lowerBound, upperBound);
74 my ivar = 0;
75 } catch (MelderError) {
76 Melder_throw (U"Linear programming: constraint not added.");
77 }
78 }
79
NUMlinprog_addConstraintCoefficient(NUMlinprog me,double coefficient)80 void NUMlinprog_addConstraintCoefficient (NUMlinprog me, double coefficient) {
81 ++ my ivar;
82 my ind [my ivar] = (int) my ivar;
83 my val [my ivar] = coefficient;
84 if (my ivar == my numberOfVariables)
85 glp_set_mat_row (my linearProgram, (int) my numberOfConstraints, (int) my numberOfVariables,
86 my ind.asArgumentToFunctionThatExpectsOneBasedArray(),
87 my val.asArgumentToFunctionThatExpectsOneBasedArray());
88 }
89
NUMlinprog_run(NUMlinprog me)90 void NUMlinprog_run (NUMlinprog me) {
91 try {
92 glp_smcp parm;
93 glp_init_smcp (& parm);
94 parm. msg_lev = GLP_MSG_OFF;
95 my status = glp_simplex (my linearProgram, & parm);
96 switch (my status) {
97 case GLP_EBADB: Melder_throw (U"Unable to start the search, because the initial basis is invalid.");
98 case GLP_ESING: Melder_throw (U"Unable to start the search, because the basis matrix is singular.");
99 case GLP_ECOND: Melder_throw (U"Unable to start the search, because the basis matrix is ill-conditioned.");
100 case GLP_EBOUND: Melder_throw (U"Unable to start the search, because some variables have incorrect bounds.");
101 case GLP_EFAIL: Melder_throw (U"Search prematurely terminated due to solver failure.");
102 case GLP_EOBJLL: Melder_throw (U"Search prematurely terminated: lower limit reached.");
103 case GLP_EOBJUL: Melder_throw (U"Search prematurely terminated: upper limit reached.");
104 case GLP_EITLIM: Melder_throw (U"Search prematurely terminated: iteration limit exceeded.");
105 case GLP_ETMLIM: Melder_throw (U"Search prematurely terminated: time limit exceeded.");
106 case GLP_ENOPFS: Melder_throw (U"The problem has no primal feasible solution.");
107 case GLP_ENODFS: Melder_throw (U"The problem has no dual feasible solution.");
108 default: break;
109 }
110 my status = glp_get_status (my linearProgram);
111 switch (my status) {
112 case GLP_INFEAS: Melder_throw (U"Solution is infeasible.");
113 case GLP_NOFEAS: Melder_throw (U"Problem has no feasible solution.");
114 case GLP_UNBND: Melder_throw (U"Problem has unbounded solution.");
115 case GLP_UNDEF: Melder_throw (U"Solution is undefined.");
116 default: break;
117 }
118 if (my status == GLP_FEAS)
119 Melder_warning (U"Linear programming solution is feasible but not optimal.");
120 } catch (MelderError) {
121 Melder_throw (U"Linear programming: not run.");
122 }
123 }
124
NUMlinprog_getPrimalValue(NUMlinprog me,integer ivar)125 double NUMlinprog_getPrimalValue (NUMlinprog me, integer ivar) {
126 return glp_get_col_prim (my linearProgram, (int) ivar);
127 }
128
129 /* End of file NUMlinprog.cpp */
130