this repo has no description
at develop 6.3 kB view raw
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 <minizinc/algorithms/min_cut.h> 23#include <minizinc/solvers/MIP/MIP_solverinstance.hh> 24 25#include <chrono> 26#include <fstream> 27#include <iomanip> 28#include <iostream> 29#include <memory> 30#include <string> 31 32using namespace std; 33 34std::string MIPWrapper::getMznLib() { return "-Glinear"; } 35 36namespace MiniZinc { 37namespace SCIPConstraints { 38 39bool check_ann_user_cut(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 check_ann_lazy_constraint(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 get_mask_cons_type(const Call* call) { 56 int mask = 0; 57 const bool fUC = check_ann_user_cut(call); 58 const bool fLC = check_ann_lazy_constraint(call); 59 if (fUC) { 60 mask |= MIPWrapper::MaskConsType_Usercut; 61 } 62 if (fLC) { 63 mask |= MIPWrapper::MaskConsType_Lazy; 64 } 65 if (!fUC && !fLC) { 66 mask |= MIPWrapper::MaskConsType_Normal; 67 } 68 return mask; 69 // return MIPWrapper::MaskConsType_Normal; // recognition fails 70} 71} // namespace SCIPConstraints 72} // namespace MiniZinc 73 74using namespace MiniZinc; 75 76void XBZCutGen::generate(const MIPWrapper::Output& slvOut, MIPWrapper::CutInput& cutsIn) { 77 assert(_pMIP); 78 const int n = static_cast<int>(varX.size()); 79 assert(n == varB.size()); 80 MIPWrapper::CutDef cut(MIPWrapper::GQ, MIPWrapper::MaskConsType_Usercut); 81 cut.addVar(varZ, -1.0); 82 for (int i = 0; i < n; ++i) { 83 const int ix = varX[i]; 84 const int ib = varB[i]; 85 assert(ix >= 0 && ix < slvOut.nCols); 86 assert(ib >= 0 && ib < slvOut.nCols); 87 const double theXi = slvOut.x[ix]; 88 const double theBi = slvOut.x[ib]; 89 const double LBXi = _pMIP->colLB[ix]; 90 const double UBXi = _pMIP->colUB[ix]; // tighter bounds from presolve? TODO 91 bool fi = (theXi + LBXi * (theBi - 1.0) - UBXi * theBi < 0.0); 92 if (fi) { 93 cut.addVar(ix, 1.0); 94 cut.addVar(ib, LBXi); 95 cut.rhs += LBXi; 96 } else { 97 cut.addVar(ib, UBXi); 98 } 99 } 100 double dViol = cut.computeViol(slvOut.x, slvOut.nCols); 101 if (dViol > 0.01) { // ?? PARAM? TODO 102 cutsIn.push_back(cut); 103 cerr << " vi" << dViol << flush; 104 // cout << cut.rmatind.size() << ' ' 105 // << cut.rhs << " cutlen, rhs. (Sense fixed to GQ) " << endl; 106 // for ( int i=0; i<cut.rmatind.size(); ++i ) 107 // cout << cut.rmatind[i] << ' '; 108 // cout << endl; 109 // for ( int i=0; i<cut.rmatind.size(); ++i ) 110 // cout << cut.rmatval[i] << ' '; 111 // cout << endl; 112 } 113} 114 115void XBZCutGen::print(ostream& os) { 116 os << varZ << '\n' << varX.size() << '\n'; 117 for (int i : varX) { 118 os << i << ' '; 119 } 120 os << endl; 121 for (int i : varB) { 122 os << i << ' '; 123 } 124 os << endl; 125} 126 127std::string SECCutGen::validate() const { 128 std::ostringstream oss; 129 /// Check that diagonal flows are 0 130 for (int i = 0; i < nN; ++i) { 131 if (_pMIP->colUB[varXij[i * nN + i]] > 0.0) { 132 oss << "SECutGen with " << nN << " cities: diagonal flow " << (i + 1) 133 << " has UB=" << _pMIP->colUB[varXij[i * nN + i]] << "\n"; 134 } 135 } 136 return oss.str(); 137} 138 139void SECCutGen::generate(const MIPWrapper::Output& slvOut, MIPWrapper::CutInput& cutsIn) { 140 assert(_pMIP); 141 /// Extract graph, converting to undirected 142 typedef map<pair<int, int>, double> TMapFlow; 143 TMapFlow mapFlow; 144 for (int i = 0; i < nN; ++i) { 145 for (int j = 0; j < nN; ++j) { 146 const double xij = slvOut.x[varXij[nN * i + j]]; 147 if (i == j) { 148 MZN_ASSERT_HARD_MSG(1e-4 > fabs(xij), 149 "circuit: X[" << (i + 1) << ", " << (j + 1) << "]==" << xij); 150 } 151 MZN_ASSERT_HARD_MSG( 152 -1e-4 < xij && 1.0 + 1e-4 > xij, // adjusted from 1e-6 to 1e-4 for CBC. 7.8.19 153 "circuit: X[" << (i + 1) << ", " << (j + 1) << "]==" << xij); 154 if (1e-4 <= xij) { 155 mapFlow[make_pair(min(i, j), max(i, j))] += xij; 156 } 157 } 158 } 159 /// Invoking Min Cut 160 // cerr << " MIN CUT... " << flush; 161 Algorithms::MinCut mc; 162 mc.nNodes = nN; 163 mc.edges.reserve(mapFlow.size()); 164 mc.weights.reserve(mapFlow.size()); 165 for (const auto& mf : mapFlow) { 166 mc.edges.push_back(mf.first); 167 mc.weights.push_back(mf.second); 168 } 169 Algorithms::MinCut::solve(); 170 /// Check if violation 171 if (mc.wMinCut <= 1.999) { 172 MIPWrapper::CutDef cut(MIPWrapper::GQ, 173 MIPWrapper::MaskConsType_Lazy | MIPWrapper::MaskConsType_Usercut); 174 cut.rhs = 1.0; 175 int nCutSize = 0; 176 constexpr int nElemPrint = 20; 177 // cerr << " CUT: [ "; 178 for (int i = 0; i < nN; ++i) { 179 if (mc.parities[i]) { 180 ++nCutSize; 181 // if ( nCutSize<=nElemPrint ) 182 // cerr << (i+1) << ", "; 183 // else if ( nCutSize==nElemPrint+1 ) 184 // cerr << "..."; 185 for (int j = 0; j < nN; ++j) { 186 if (!mc.parities[j]) { 187 cut.addVar(varXij[nN * i + j], 1.0); 188 } 189 } 190 } 191 } 192 // cerr << "]. " << flush; 193 double dViol = cut.computeViol(slvOut.x, slvOut.nCols); 194 if (dViol > 0.0001) { // ?? PARAM? TODO. See also min cut value required 195 cutsIn.push_back(cut); 196 /* cerr << " SEC: viol=" << dViol 197 << " N NODES: " << nN 198 << " |X|: : " << nCutSize 199 << flush; */ 200 } else { 201 MZN_ASSERT_HARD_MSG(0, " SEC cut: N nodes = " << nN << ": violation = " << dViol 202 << ": too small compared to the min-cut value " 203 << (2.0 - mc.wMinCut)); 204 } 205 } 206} 207 208void SECCutGen::print(ostream& /*os*/) {}