this repo has no description
dekker.one
1// * -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */
2
3/*
4 * Main authors:
5 * Gleb Belov <gleb.belov@monash.edu>
6 */
7
8/* This Source Code Form is subject to the terms of the Mozilla Public
9 * License, v. 2.0. If a copy of the MPL was ! distributed with this
10 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
11
12/* This is the main file for a mzn-cplex solver using a unified
13 * linearization module && a flexible flattener-to-solver interface
14 */
15
16/// TODO Quadratic terms, even CBC
17
18#ifdef _MSC_VER
19#define _CRT_SECURE_NO_WARNINGS
20#endif
21
22#include <chrono>
23#include <fstream>
24#include <iomanip>
25#include <iostream>
26#include <memory>
27#include <string>
28
29using namespace std;
30
31#include <minizinc/algorithms/min_cut.h>
32#include <minizinc/solvers/MIP/MIP_solverinstance.hh>
33
34std::string MIP_wrapper::getMznLib() { return "-Glinear"; }
35
36namespace MiniZinc {
37namespace SCIPConstraints {
38
39bool CheckAnnUserCut(const Call* call) {
40 if (!call->ann().isEmpty()) {
41 if (call->ann().contains(constants().ann.user_cut)) {
42 return true;
43 }
44 }
45 return false;
46}
47bool CheckAnnLazyConstraint(const Call* call) {
48 if (!call->ann().isEmpty()) {
49 if (call->ann().contains(constants().ann.lazy_constraint)) {
50 return true;
51 }
52 }
53 return false;
54}
55int GetMaskConsType(const Call* call) {
56 int mask = 0;
57 const bool fUC = CheckAnnUserCut(call);
58 const bool fLC = CheckAnnLazyConstraint(call);
59 if (fUC) {
60 mask |= MIP_wrapper::MaskConsType_Usercut;
61 }
62 if (fLC) {
63 mask |= MIP_wrapper::MaskConsType_Lazy;
64 }
65 if (!fUC && !fLC) mask |= MIP_wrapper::MaskConsType_Normal;
66 return mask;
67 // return MIP_wrapper::MaskConsType_Normal; // recognition fails
68}
69} // namespace SCIPConstraints
70} // namespace MiniZinc
71
72using namespace MiniZinc;
73
74void XBZCutGen::generate(const MIP_wrapper::Output& slvOut, MIP_wrapper::CutInput& cutsIn) {
75 assert(pMIP);
76 const int n = static_cast<int>(varX.size());
77 assert(n == varB.size());
78 MIP_wrapper::CutDef cut(MIP_wrapper::GQ, MIP_wrapper::MaskConsType_Usercut);
79 cut.addVar(varZ, -1.0);
80 for (int i = 0; i < n; ++i) {
81 const int ix = varX[i];
82 const int ib = varB[i];
83 assert(ix >= 0 && ix < slvOut.nCols);
84 assert(ib >= 0 && ib < slvOut.nCols);
85 const double theXi = slvOut.x[ix];
86 const double theBi = slvOut.x[ib];
87 const double LBXi = pMIP->colLB[ix];
88 const double UBXi = pMIP->colUB[ix]; // tighter bounds from presolve? TODO
89 bool fi = (theXi + LBXi * (theBi - 1.0) - UBXi * theBi < 0.0);
90 if (fi) {
91 cut.addVar(ix, 1.0);
92 cut.addVar(ib, LBXi);
93 cut.rhs += LBXi;
94 } else {
95 cut.addVar(ib, UBXi);
96 }
97 }
98 double dViol = cut.computeViol(slvOut.x, slvOut.nCols);
99 if (dViol > 0.01) { // ?? PARAM? TODO
100 cutsIn.push_back(cut);
101 cerr << " vi" << dViol << flush;
102 // cout << cut.rmatind.size() << ' '
103 // << cut.rhs << " cutlen, rhs. (Sense fixed to GQ) " << endl;
104 // for ( int i=0; i<cut.rmatind.size(); ++i )
105 // cout << cut.rmatind[i] << ' ';
106 // cout << endl;
107 // for ( int i=0; i<cut.rmatind.size(); ++i )
108 // cout << cut.rmatval[i] << ' ';
109 // cout << endl;
110 }
111}
112
113void XBZCutGen::print(ostream& os) {
114 os << varZ << '\n' << varX.size() << '\n';
115 for (int i = 0; i < varX.size(); ++i) os << varX[i] << ' ';
116 os << endl;
117 for (int i = 0; i < varB.size(); ++i) os << varB[i] << ' ';
118 os << endl;
119}
120
121std::string SECCutGen::validate() const {
122 std::ostringstream oss;
123 /// Check that diagonal flows are 0
124 for (int i = 0; i < nN; ++i)
125 if (pMIP->colUB[varXij[i * nN + i]] > 0.0)
126 oss << "SECutGen with " << nN << " cities: diagonal flow " << (i + 1)
127 << " has UB=" << pMIP->colUB[varXij[i * nN + i]] << "\n";
128 return oss.str();
129}
130
131void SECCutGen::generate(const MIP_wrapper::Output& slvOut, MIP_wrapper::CutInput& cutsIn) {
132 assert(pMIP);
133 /// Extract graph, converting to undirected
134 typedef map<pair<int, int>, double> TMapFlow;
135 TMapFlow mapFlow;
136 for (int i = 0; i < nN; ++i) {
137 for (int j = 0; j < nN; ++j) {
138 const double xij = slvOut.x[varXij[nN * i + j]];
139 if (i == j)
140 MZN_ASSERT_HARD_MSG(1e-4 > fabs(xij),
141 "circuit: X[" << (i + 1) << ", " << (j + 1) << "]==" << xij);
142 MZN_ASSERT_HARD_MSG(
143 -1e-4 < xij && 1.0 + 1e-4 > xij, // adjusted from 1e-6 to 1e-4 for CBC. 7.8.19
144 "circuit: X[" << (i + 1) << ", " << (j + 1) << "]==" << xij);
145 if (1e-4 <= xij) {
146 mapFlow[make_pair(min(i, j), max(i, j))] += xij;
147 }
148 }
149 }
150 /// Invoking Min Cut
151 // cerr << " MIN CUT... " << flush;
152 Algorithms::MinCut mc;
153 mc.nNodes = nN;
154 mc.edges.reserve(mapFlow.size());
155 mc.weights.reserve(mapFlow.size());
156 for (const auto& mf : mapFlow) {
157 mc.edges.push_back(mf.first);
158 mc.weights.push_back(mf.second);
159 }
160 mc.solve();
161 /// Check if violation
162 if (mc.wMinCut <= 1.999) {
163 MIP_wrapper::CutDef cut(MIP_wrapper::GQ,
164 MIP_wrapper::MaskConsType_Lazy | MIP_wrapper::MaskConsType_Usercut);
165 cut.rhs = 1.0;
166 int nCutSize = 0;
167 constexpr int nElemPrint = 20;
168 // cerr << " CUT: [ ";
169 for (int i = 0; i < nN; ++i)
170 if (mc.parities[i]) {
171 ++nCutSize;
172 // if ( nCutSize<=nElemPrint )
173 // cerr << (i+1) << ", ";
174 // else if ( nCutSize==nElemPrint+1 )
175 // cerr << "...";
176 for (int j = 0; j < nN; ++j)
177 if (!mc.parities[j]) {
178 cut.addVar(varXij[nN * i + j], 1.0);
179 }
180 }
181 // cerr << "]. " << flush;
182 double dViol = cut.computeViol(slvOut.x, slvOut.nCols);
183 if (dViol > 0.0001) { // ?? PARAM? TODO. See also min cut value required
184 cutsIn.push_back(cut);
185 /* cerr << " SEC: viol=" << dViol
186 << " N NODES: " << nN
187 << " |X|: : " << nCutSize
188 << flush; */
189 } else {
190 MZN_ASSERT_HARD_MSG(0, " SEC cut: N nodes = " << nN << ": violation = " << dViol
191 << ": too small compared to the min-cut value "
192 << (2.0 - mc.wMinCut));
193 }
194 }
195}
196
197void SECCutGen::print(ostream&) {}