this repo has no description
dekker.one
1/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */
2
3/* This Source Code Form is subject to the terms of the Mozilla Public
4 * License, v. 2.0. If a copy of the MPL was not distributed with this
5 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
6
7#include <minizinc/solvers/nl/nl_components.hh>
8#include <minizinc/solvers/nl/nl_file.hh>
9
10#include <utility>
11
12using namespace std;
13
14namespace MiniZinc {
15
16/* *** *** *** NLBound *** *** *** */
17
18// Constructors
19NLBound::NLBound(Bound tag, double lb, double ub) : tag(tag), lb(lb), ub(ub) {}
20
21NLBound NLBound::makeBounded(double lb, double ub) {
22 if (lb == ub) {
23 return makeEqual(lb);
24 }
25 assert(lb < ub);
26 return NLBound(LB_UB, lb, ub);
27}
28
29NLBound NLBound::makeUBBounded(double ub) { return NLBound(UB, 0, ub); }
30
31NLBound NLBound::makeLBBounded(double lb) { return NLBound(LB, lb, 0); }
32
33NLBound NLBound::makeNoBound() { return NLBound(NONE, 0, 0); }
34
35NLBound NLBound::makeEqual(double val) { return NLBound(EQ, val, val); }
36
37/** Update the lower bound */
38void NLBound::updateLB(double new_lb) {
39 switch (tag) {
40 // LB <= var <= UB. Same tag
41 case NLBound::LB_UB: {
42 assert(new_lb <= ub);
43 if (new_lb > lb) {
44 lb = new_lb;
45 }
46 break;
47 }
48 // var <= UB. Update tag
49 case NLBound::UB: {
50 assert(new_lb <= ub);
51 tag = LB_UB;
52 lb = new_lb;
53 break;
54 }
55 // LB <= var. Same tag
56 case NLBound::LB: {
57 if (new_lb > lb) {
58 lb = new_lb;
59 }
60 break;
61 }
62 // No bound. Update tag
63 case NLBound::NONE: {
64 tag = LB;
65 lb = new_lb;
66 break;
67 }
68 // LB = var = UB. Should not happen
69 case NLBound::EQ: {
70 should_not_happen(
71 "Updating a bound already set to \"equals\". We only allow tightening update.");
72 }
73 }
74}
75
76/** Update the upper bound */
77void NLBound::updateUB(double new_ub) {
78 switch (tag) {
79 // LB <= var <= UB. Same tag
80 case NLBound::LB_UB: {
81 assert(lb <= new_ub);
82 if (new_ub < ub) {
83 ub = new_ub;
84 }
85 break;
86 }
87 // var <= UB. Same tag
88 case NLBound::UB: {
89 if (new_ub < ub) {
90 ub = new_ub;
91 }
92 break;
93 }
94 // LB <= var. Update tag
95 case NLBound::LB: {
96 assert(lb <= new_ub);
97 tag = LB_UB;
98 ub = new_ub;
99 break;
100 }
101 // No bound. Update tag
102 case NLBound::NONE: {
103 tag = UB;
104 ub = new_ub;
105 break;
106 }
107 // LB = var = UB. Should not happen
108 case NLBound::EQ: {
109 should_not_happen(
110 "Updating a bound already set to \"equals\". We only allow tightening update.");
111 }
112 }
113}
114
115/** Update equal bound */
116void NLBound::updateEq(double new_eq) {
117 switch (tag) {
118 // LB <= var <= UB. Update tag
119 case NLBound::LB_UB: {
120 assert(lb <= new_eq && new_eq <= ub);
121 tag = EQ;
122 lb = new_eq;
123 ub = new_eq;
124 }
125 // var <= UB. Update tag
126 case NLBound::UB: {
127 assert(new_eq <= ub);
128 tag = EQ;
129 lb = new_eq;
130 ub = new_eq;
131 }
132 // LB <= var. Update tag
133 case NLBound::LB: {
134 assert(lb <= new_eq);
135 tag = EQ;
136 lb = new_eq;
137 ub = new_eq;
138 }
139 // No bound. Update tag
140 case NLBound::NONE: {
141 tag = EQ;
142 lb = new_eq;
143 ub = new_eq;
144 }
145 // LB = var = UB. Can happen only if we "update" with the same bound.
146 case NLBound::EQ: {
147 assert(lb == new_eq);
148 }
149 }
150}
151
152/** Printing with a name. */
153ostream& NLBound::printToStream(ostream& os, const string& vname) const {
154 switch (tag) {
155 case LB_UB: {
156 os << "0 " << lb << " " << ub << " # " << lb << " =< " << vname << " =< " << ub;
157 break;
158 }
159 case UB: {
160 os << "1 " << ub << " # " << vname << " =< " << ub;
161 break;
162 }
163 case LB: {
164 os << "2 " << lb << " # " << lb << " =< " << vname;
165 break;
166 }
167 case NONE: {
168 os << "3"
169 << " # No constraint";
170 break;
171 }
172 case EQ: {
173 os << "4 " << lb << " # " << vname << " = " << lb;
174 break;
175 }
176 }
177 return os;
178}
179
180/** Default printing */
181ostream& NLBound::printToStream(ostream& os) const {
182 printToStream(os, "body");
183 return os;
184}
185
186/* *** *** *** NLVar *** *** *** */
187
188NLVar NLVar::copyWithBound(NLBound bound) const {
189 NLVar v = NLVar(*this); // copy constructor
190 v.bound = bound;
191 return v;
192}
193
194/* *** *** *** NLArray *** *** *** */
195
196/* *** *** *** NLToken *** *** *** */
197
198const char* NLToken::getName(OpCode oc) {
199 switch (oc) {
200 case OpCode::OPPLUS:
201 return "OPPLUS";
202 case OpCode::OPMINUS:
203 return "OPMINUS";
204 case OpCode::OPMULT:
205 return "OPMULT";
206 case OpCode::OPDIV:
207 return "OPDIV";
208 case OpCode::OPREM:
209 return "OPREM";
210 case OpCode::OPPOW:
211 return "OPPOW";
212 case OpCode::OPLESS:
213 return "OPLESS";
214 case OpCode::FLOOR:
215 return "FLOOR";
216 case OpCode::CEIL:
217 return "CEIL";
218 case OpCode::ABS:
219 return "ABS";
220 case OpCode::OPUMINUS:
221 return "OPUMINUS";
222 case OpCode::OPOR:
223 return "OPOR";
224 case OpCode::OPAND:
225 return "OPAND";
226 case OpCode::LT:
227 return "LT";
228 case OpCode::LE:
229 return "LE";
230 case OpCode::EQ:
231 return "EQ";
232 case OpCode::GE:
233 return "GE";
234 case OpCode::GT:
235 return "GT";
236 case OpCode::NE:
237 return "NE";
238 case OpCode::OPNOT:
239 return "OPNOT";
240 case OpCode::OPIFnl:
241 return "OPIFnl";
242 case OpCode::OP_tanh:
243 return "OP_tanh";
244 case OpCode::OP_tan:
245 return "OP_tan";
246 case OpCode::OP_sqrt:
247 return "OP_sqrt";
248 case OpCode::OP_sinh:
249 return "OP_sinh";
250 case OpCode::OP_sin:
251 return "OP_sin";
252 case OpCode::OP_log10:
253 return "OP_log10";
254 case OpCode::OP_log:
255 return "OP_log";
256 case OpCode::OP_exp:
257 return "OP_exp";
258 case OpCode::OP_cosh:
259 return "OP_cosh";
260 case OpCode::OP_cos:
261 return "OP_cos";
262 case OpCode::OP_atanh:
263 return "OP_atanh";
264 case OpCode::OP_atan2:
265 return "OP_atan2";
266 case OpCode::OP_atan:
267 return "OP_atan";
268 case OpCode::OP_asinh:
269 return "OP_asinh";
270 case OpCode::OP_asin:
271 return "OP_asin";
272 case OpCode::OP_acosh:
273 return "OP_acosh";
274 case OpCode::OP_acos:
275 return "OP_acos";
276 case OpCode::OPintDIV:
277 return "OPintDIV";
278 case OpCode::OPprecision:
279 return "OPprecision";
280 case OpCode::OPround:
281 return "OPround";
282 case OpCode::OPtrunc:
283 return "OPtrunc";
284 case OpCode::OPATLEAST:
285 return "OPATLEAST";
286 case OpCode::OPATMOST:
287 return "OPATMOST";
288 case OpCode::OPPLTERM:
289 return "OPPLTERM";
290 case OpCode::OPIFSYM:
291 return "OPIFSYM";
292 case OpCode::OPEXACTLY:
293 return "OPEXACTLY";
294 case OpCode::OPNOTATLEAST:
295 return "OPNOTATLEAST";
296 case OpCode::OPNOTATMOST:
297 return "OPNOTATMOST";
298 case OpCode::OPNOTEXACTLY:
299 return "OPNOTEXACTLY";
300 case OpCode::OPIMPELSE:
301 return "OPIMPELSE";
302 case OpCode::OP_IFF:
303 return "OP_IFF";
304 case OpCode::OPSOMESAME:
305 return "OPSOMESAME";
306 case OpCode::OP1POW:
307 return "OP1POW";
308 case OpCode::OP2POW:
309 return "OP2POW";
310 case OpCode::OPCPOW:
311 return "OPCPOW";
312 case OpCode::OPFUNCALL:
313 return "OPFUNCALL";
314 case OpCode::OPNUM:
315 return "OPNUM";
316 case OpCode::OPHOL:
317 return "OPHOL";
318 case OpCode::OPVARVAL:
319 return "OPVARVAL";
320 case OpCode::N_OPS:
321 return "N_OPS";
322 default:
323 assert(false);
324 return nullptr;
325 }
326};
327
328const char* NLToken::getName(MOpCode moc) {
329 switch (moc) {
330 case MOpCode::MINLIST:
331 return "MINLIST";
332 case MOpCode::MAXLIST:
333 return "MAXLIST";
334 case MOpCode::OPSUMLIST:
335 return "OPSUMLIST";
336 case MOpCode::OPCOUNT:
337 return "OPCOUNT";
338 case MOpCode::OPNUMBEROF:
339 return "OPNUMBEROF";
340 case MOpCode::OPNUMBEROFs:
341 return "OPNUMBEROFs";
342 case MOpCode::ANDLIST:
343 return "ANDLIST";
344 case MOpCode::ORLIST:
345 return "ORLIST";
346 case MOpCode::OPALLDIFF:
347 return "OPALLDIFF";
348 default:
349 assert(false);
350 return nullptr;
351 }
352}
353
354NLToken NLToken::n(double value) {
355 NLToken tok;
356 tok.kind = Kind::NUMERIC;
357 tok.numericValue = value;
358 return tok;
359}
360
361NLToken NLToken::v(string vname) {
362 NLToken tok;
363 tok.kind = Kind::VARIABLE;
364 tok.str = std::move(vname);
365 return tok;
366}
367
368NLToken NLToken::o(OpCode opc) {
369 NLToken tok;
370 tok.kind = Kind::OP;
371 tok.oc = opc;
372 return tok;
373}
374
375NLToken NLToken::mo(MOpCode mopc, int nb) {
376 NLToken tok;
377 tok.kind = Kind::MOP;
378 tok.moc = mopc;
379 tok.argCount = nb;
380 return tok;
381}
382
383bool NLToken::isVariable() const { return kind == VARIABLE; }
384
385bool NLToken::isConstant() const { return kind == NUMERIC; }
386
387ostream& NLToken::printToStream(ostream& os, const NLFile& nl_file) const {
388 switch (kind) {
389 case Kind::NUMERIC: {
390 os << "n" << numericValue;
391 break;
392 }
393
394 case Kind::VARIABLE: {
395 os << "v" << nl_file.variableIndexes.at(str) << " # " << str;
396 break;
397 }
398
399 case Kind::STRING: {
400 should_not_happen("NL string token (Kind::STRING) not implemented");
401 }
402
403 case Kind::FUNCALL: {
404 should_not_happen("NL function call token (Kind::FUNCALL) not implemented");
405 }
406
407 case Kind::OP: {
408 os << "o" << oc << " # " << getName(oc);
409 break;
410 }
411
412 case Kind::MOP: {
413 os << "o" << moc << " # " << getName(moc) << " " << argCount << endl;
414 os << argCount;
415 break;
416 }
417
418 default:
419 should_not_happen("Unknown token kind: " << kind);
420 }
421
422 return os;
423}
424
425/* *** *** *** NLAlgCons *** *** *** */
426
427/** Method to build the var_coeff vector. */
428void NLAlgCons::setJacobian(const vector<string>& vnames, const vector<double>& coeffs,
429 NLFile* nl_file) {
430 assert(vnames.size() == coeffs.size());
431 for (int i = 0; i < vnames.size(); ++i) {
432 string vn = vnames[i];
433 nl_file->variables.at(vn).jacobianCount++;
434 jacobian.emplace_back(vn, coeffs[i]);
435 }
436}
437
438/** A constraint is considered linear if the expressionGraph is empty. */
439bool NLAlgCons::isLinear() const { return expressionGraph.empty(); }
440
441/** Printing. */
442ostream& NLAlgCons::printToStream(ostream& os, const NLFile& nl_file) const {
443 int idx = nl_file.constraintIndexes.at(name);
444
445 // Print the 'C' segment: if no expression graph, print "n0".
446 os << "C" << idx << " # Non linear part of " << name << endl;
447 if (expressionGraph.empty()) {
448 os << "n0 # No non linear part coded as the value '0'" << endl;
449 } else {
450 for (const auto& t : expressionGraph) {
451 t.printToStream(os, nl_file);
452 os << endl;
453 }
454 }
455
456 // Print the 'J' segment if present.
457 if (!jacobian.empty()) {
458 os << "J" << idx << " " << jacobian.size() << " # Linear part of " << name << endl;
459 for (const auto& vn_coef : jacobian) {
460 os << nl_file.variableIndexes.at(vn_coef.first) << " " << vn_coef.second << " # "
461 << vn_coef.first << endl;
462 }
463 }
464
465 return os;
466}
467
468/* *** *** *** NLLogicalCons *** *** *** */
469
470/** Printing. */
471ostream& NLLogicalCons::printToStream(ostream& os, const NLFile& nl_file) const {
472 os << "L" << index << " # Logical constraint " << name << endl;
473 for (const auto& t : expressionGraph) {
474 t.printToStream(os, nl_file);
475 os << endl;
476 }
477
478 return os;
479}
480
481/* *** *** *** NLHeader *** *** *** */
482
483/** Printing the header.
484 * The header is composed of then lines that we describe as we proceed.
485 * A '#' starts a comment until the end of the line. However, it cannot be a line on its own!*/
486ostream& NLHeader::printToStream(ostream& os, const NLFile& nl_file) {
487 // 1st line:
488 // 'g': file will be in text format
489 // other numbers: as given in the doc (no other explanation...)
490 os << "g3 1 1 0" << endl;
491
492 // 2nd line:
493 os << nl_file.variables.size() << " " // Total number of variables
494 << nl_file.constraints.size()
495 << " " // Total number of algebraic constraint (including 'range' and 'eq')
496 << 1 << " " // Always 1 objective
497 << nl_file.algConsRangeCount << " " // Number of algebraic range constraints
498 << nl_file.algConsEqCount << " " // Number of algebraic eq constraints
499 << nl_file.logicalConstraints.size() << " " // Number of logical constraints
500 << "# Total nb of: variables, algebraic constraints, objectives, ranges, eqs, logical "
501 "constraints"
502 << endl;
503
504 // 3rd line: Nonlinear and complementary information
505 os << nl_file.cnames_nl_general.size() << " " // Non linear constraints
506 << (nl_file.objective.isLinear() ? 0 : 1) << " " // Non linear objective
507 << "# Nb of nonlinear constraints, nonlinar objectives." << endl;
508 /* This was found in the online source of the ASL parser, but is not produce in our ampl tests.
509 * If needed, should be put on the same line
510 << nb_complementarity_linear_conditions << " "
511 << nb_complementarity_nonlinear_conditions << " "
512 << nb_complementarity_items_double_inequalities << " "
513 << nb_complemented_vars_non_zero_lowerbound << " "
514 << "Nb of complementary: linear & nonlinear conditions, double inequalities, vars with non-0 lower
515 bound."
516 << endl;
517 */
518
519 // 4th line: Network constraints
520 os << nl_file.cnames_nl_network.size() << " " // Number of nonlinear network constraints
521 << nl_file.cnames_lin_network.size() << " " // Number of linear network constraints
522 << "# Nb of network constraints: nonlinear, linear." << endl;
523
524 // 5th line: nonlinear variables:
525 os << nl_file.lvcCount() << " " // Nb of nonlinear vars in constraints
526 << nl_file.lvoCount() << " " // Nb of nonlinear vars in objectives
527 << nl_file.lvbCount() << " " // Nb of nonlinear vars in both
528 << "# Nb of non linear vars in: constraints, objectives, both." << endl;
529
530 // 6th line:
531 os << nl_file.wvCount() << " " // Nb of linear network vars
532 << "0"
533 << " " // Nb of functions. Not Implemented
534 << "0 1 "
535 << "# Nb of: linear network vars, functions. Floating point arithmetic mode (TEXT == 0). "
536 "Flag: if 1, add .sol suffixe."
537 << endl;
538
539 // 7th line: discrete variables
540 os << nl_file.bvCount() << " " // Nb of linear binary vars
541 << nl_file.ivCount() << " " // Nb of linear integer vars
542 << nl_file.lvbiCount() << " " // Nb of nonlinear integer vars in both
543 << nl_file.lvciCount() << " " // Nb of nonlinear integer vars in constraints only
544 << nl_file.lvoiCount() << " " // Nb of nonlinear integer vars in objectives only
545 << "# Nb of linear vars: binary, integer (non binary). "
546 << "Nb of nonlinear integer vars in: both, constraints only, objectives only." << endl;
547
548 // 8th line: non zeros
549 os << nl_file.jacobianCount() << " " // Nb of nonzero in jacobian
550 << nl_file.objective.gradientCount() << " " // Nb of nonzero in gradient
551 << "# Nb of non zeros in: jacobian, objective gradients." << endl;
552
553 // 9th line: name length. Our tests always produce 0...
554 os << "0"
555 << " " // Max constraint name length (??)
556 << "0"
557 << " " // Max var name length (??)
558 << "# Longest name among: contraints' name, vars' name." << endl;
559
560 // 10th line: common expressions. Not enough infor for now...
561 os << "0"
562 << " " // Nb common exprs in both
563 << "0"
564 << " " // Nb common exprs in constraints only
565 << "0"
566 << " " // Nb common exprs in objectives only
567 << "0"
568 << " " // Nb common exprs in single constraint only
569 << "0"
570 << " " // Nb common exprs in single objectives only
571 << "# Nb of common expressions in: both, constraints only, objectives only, single "
572 "constraint, single objective.";
573
574 return os;
575}
576
577/* *** *** *** NLObjective *** *** *** */
578
579/** Gradient count. */
580int NLObjective::gradientCount() const { return gradient.size(); }
581
582/** Set the gradient. */
583void NLObjective::setGradient(const vector<string>& vnames, const vector<double>& coeffs) {
584 assert(vnames.size() == coeffs.size());
585 for (int i = 0; i < vnames.size(); ++i) {
586 string vn = vnames[i];
587 gradient.emplace_back(vn, coeffs[i]);
588 }
589}
590
591/** A objective is considered as linear if its expression graph is non empty. */
592bool NLObjective::isDefined() const { return minmax != UNDEF; }
593
594bool NLObjective::isLinear() const { return expressionGraph.empty(); }
595
596bool NLObjective::isOptimisation() const { return minmax >= MINIMIZE; }
597
598/** Printing. */
599ostream& NLObjective::printToStream(ostream& os, const NLFile& nl_file) const {
600 if (minmax != UNDEF) {
601 if (minmax == SATISFY) {
602 os << "O0 0 # Satisfy objectif implemented as 'minimize 0'" << endl;
603 os << "n0" << endl;
604 } else {
605 os << "O0 " << minmax << " # Objectif (0: minimize, 1: maximize)" << endl;
606 if (expressionGraph.empty()) {
607 os << "n0 # No expression graph" << endl;
608 } else {
609 for (const auto& tok : expressionGraph) {
610 tok.printToStream(os, nl_file) << endl;
611 }
612 }
613 // Print gradient
614 if (!gradient.empty()) {
615 os << "G0 " << gradient.size() << " # Objective Linear part" << endl;
616 for (const auto& vn_coef : gradient) {
617 os << nl_file.variableIndexes.at(vn_coef.first) << " " << vn_coef.second << " # "
618 << vn_coef.first << endl;
619 }
620 }
621 }
622 }
623
624 return os;
625}
626} // namespace MiniZinc