|
| 1 | +// Copyright (C) 2022 Mathieu Dutour Sikiric <[email protected]> |
| 2 | +#include "MAT_Matrix.h" |
| 3 | +int main(int argc, char *argv[]) { |
| 4 | + try { |
| 5 | + if (argc != 2) { |
| 6 | + fprintf(stderr, "Number of argument is = %d\n", argc); |
| 7 | + fprintf(stderr, "ComputeChromaticLowerBounds [inputAdjMat]\n"); |
| 8 | + return -1; |
| 9 | + } |
| 10 | + // |
| 11 | + std::string eFile = argv[1]; |
| 12 | + auto GetListDegree = [&]() -> std::vector<int> { |
| 13 | + std::ifstream is(eFile); |
| 14 | + int nbVert, eAdj; |
| 15 | + is >> nbVert; |
| 16 | + std::vector<int> ListDeg(nbVert); |
| 17 | + for (int iVert = 0; iVert < nbVert; iVert++) { |
| 18 | + int nbAdj; |
| 19 | + is >> nbAdj; |
| 20 | + for (int iAdj = 0; iAdj < nbAdj; iAdj++) |
| 21 | + is >> eAdj; |
| 22 | + ListDeg[iVert] = nbAdj; |
| 23 | + } |
| 24 | + return ListDeg; |
| 25 | + }; |
| 26 | + // |
| 27 | + std::vector<int> ListDeg = GetListDegree(); |
| 28 | + int nnz = 0; |
| 29 | + for (auto &eVal : ListDeg) |
| 30 | + nnz += eVal; |
| 31 | + // |
| 32 | + using T2 = Eigen::Triplet<double>; |
| 33 | + std::ifstream is(eFile); |
| 34 | + int nbVert; |
| 35 | + is >> nbVert; |
| 36 | + std::vector<T2> tripletList_A(nnz); |
| 37 | + std::vector<T2> tripletList_D(nbVert); |
| 38 | + int iNNZ = 0; |
| 39 | + double one_val = 1.0; |
| 40 | + for (int iVert = 0; iVert < nbVert; iVert++) { |
| 41 | + int nbAdj; |
| 42 | + is >> nbAdj; |
| 43 | + for (int iAdj = 0; iAdj < nbAdj; iAdj++) { |
| 44 | + int eAdj; |
| 45 | + is >> eAdj; |
| 46 | + tripletList_A[iNNZ] = T2(iVert, eAdj, one_val); |
| 47 | + iNNZ++; |
| 48 | + } |
| 49 | + tripletList_D[iVert] = T2(iVert, iVert, nbAdj); |
| 50 | + } |
| 51 | + MySparseMatrix<double> A(nbVert, nbVert); |
| 52 | + A.setFromTriplets(tripletList_A.begin(), tripletList_A.end()); |
| 53 | + std::cerr << "We have A matrix\n"; |
| 54 | + MySparseMatrix<double> D(nbVert, nbVert); |
| 55 | + D.setFromTriplets(tripletList_D.begin(), tripletList_D.end()); |
| 56 | + std::cerr << "We have D matrix\n"; |
| 57 | + // |
| 58 | + auto GetEigenvalues = |
| 59 | + [&](MySparseMatrix<double> const &SpMat) -> MyVector<double> { |
| 60 | + Eigen::SelfAdjointEigenSolver<MySparseMatrix<double>> eig(SpMat); |
| 61 | + MyVector<double> ListEig = eig.eigenvalues(); |
| 62 | + MyVector<double> ListEigRev(nbVert); |
| 63 | + for (int iVert = 0; iVert < nbVert; iVert++) |
| 64 | + ListEigRev(iVert) = ListEig(nbVert - 1 - iVert); |
| 65 | + return ListEigRev; |
| 66 | + }; |
| 67 | + // |
| 68 | + MyVector<double> mu = GetEigenvalues(A); |
| 69 | + std::cerr << "We have mu eigenvalues\n"; |
| 70 | + MyVector<double> delta = GetEigenvalues(A + D); |
| 71 | + std::cerr << "We have delta eigenvalues\n"; |
| 72 | + MyVector<double> theta = GetEigenvalues(D - A); |
| 73 | + std::cerr << "We have theta eigenvalues\n"; |
| 74 | + int n = nbVert; |
| 75 | + // |
| 76 | + double BestLowerBound = 0; |
| 77 | + std::string RealizingBound = "unset"; |
| 78 | + auto UpdateBounds = [&](double const &eBound, |
| 79 | + std::string const &RealBound) -> void { |
| 80 | + if (eBound > BestLowerBound) { |
| 81 | + BestLowerBound = eBound; |
| 82 | + RealizingBound = RealBound; |
| 83 | + } |
| 84 | + }; |
| 85 | + // |
| 86 | + double LowerBoundHoffman = 1 + mu(0) / (-mu(n - 1)); |
| 87 | + UpdateBounds(LowerBoundHoffman, "Hoffman bound"); |
| 88 | + std::cerr << "Hoffman lower bound = " << LowerBoundHoffman << "\n"; |
| 89 | + // |
| 90 | + double LowerBoundNifikorov = 1 + mu(0) / (theta(0) - mu(0)); |
| 91 | + UpdateBounds(LowerBoundNifikorov, "Nifikorov bound"); |
| 92 | + std::cerr << "Nifikorov lower bound = " << LowerBoundNifikorov << "\n"; |
| 93 | + // |
| 94 | + double LowerBoundKotolina1 = 1 + mu(0) / (mu(0) - delta(0) + theta(0)); |
| 95 | + UpdateBounds(LowerBoundKotolina1, "Kotolina1 bound"); |
| 96 | + std::cerr << "Kotolina lower bound 1 = " << LowerBoundKotolina1 << "\n"; |
| 97 | + // |
| 98 | + double LowerBoundKotolina2 = |
| 99 | + 1 + mu(0) / (mu(0) - delta(n - 1) + theta(n - 1)); |
| 100 | + UpdateBounds(LowerBoundKotolina2, "Kotolina2 bound"); |
| 101 | + std::cerr << "Kotolina lower bound 2 = " << LowerBoundKotolina2 << "\n"; |
| 102 | + // |
| 103 | + for (int m = 1; m <= n; m++) { |
| 104 | + double sum1; |
| 105 | + double sum2; |
| 106 | + // |
| 107 | + sum1 = 0; |
| 108 | + sum2 = 0; |
| 109 | + for (int i = 1; i <= m; i++) { |
| 110 | + sum1 += mu(i - 1); |
| 111 | + sum2 += mu(n + 1 - i - 1); |
| 112 | + } |
| 113 | + double UnifiedLower1 = 1 + sum1 / (-sum2); |
| 114 | + UpdateBounds(UnifiedLower1, |
| 115 | + "Unified lower bound 1: m=" + std::to_string(m)); |
| 116 | + std::cerr << "Elphick Wocjan unified lower bound 1 = " << UnifiedLower1 |
| 117 | + << "\n"; |
| 118 | + // |
| 119 | + sum1 = 0; |
| 120 | + sum2 = 0; |
| 121 | + for (int i = 1; i <= m; i++) { |
| 122 | + sum1 += mu(i - 1); |
| 123 | + sum2 += theta(i - 1) - mu(i - 1); |
| 124 | + } |
| 125 | + double UnifiedLower2 = 1 + sum1 / sum2; |
| 126 | + UpdateBounds(UnifiedLower2, |
| 127 | + "Unified lower bound 2: m=" + std::to_string(m)); |
| 128 | + std::cerr << "Elphick Wocjan unified lower bound 2 = " << UnifiedLower2 |
| 129 | + << "\n"; |
| 130 | + // |
| 131 | + sum1 = 0; |
| 132 | + sum2 = 0; |
| 133 | + for (int i = 1; i <= m; i++) { |
| 134 | + sum1 += mu(i - 1); |
| 135 | + sum2 += mu(i - 1) - delta(i - 1) + theta(i - 1); |
| 136 | + } |
| 137 | + double UnifiedLower3 = 1 + sum1 / sum2; |
| 138 | + UpdateBounds(UnifiedLower3, |
| 139 | + "Unified lower bound 3: m=" + std::to_string(m)); |
| 140 | + std::cerr << "Elphick Wocjan unified lower bound 3 = " << UnifiedLower3 |
| 141 | + << "\n"; |
| 142 | + // |
| 143 | + sum1 = 0; |
| 144 | + sum2 = 0; |
| 145 | + for (int i = 1; i <= m; i++) { |
| 146 | + sum1 += mu(i - 1); |
| 147 | + sum2 += mu(i - 1) - delta(n - i) + theta(n - i); |
| 148 | + } |
| 149 | + double UnifiedLower4 = 1 + sum1 / sum2; |
| 150 | + UpdateBounds(UnifiedLower4, |
| 151 | + "Unified lower bound 4: m=" + std::to_string(m)); |
| 152 | + std::cerr << "Elphick Wocjan unified lower bound 4 = " << UnifiedLower4 |
| 153 | + << "\n"; |
| 154 | + } |
| 155 | + // |
| 156 | + double splus = 0; |
| 157 | + double sminus = 0; |
| 158 | + for (int i = 0; i < n; i++) { |
| 159 | + if (mu(i) > 0) |
| 160 | + splus += mu(i) * mu(i); |
| 161 | + if (mu(i) < 0) |
| 162 | + sminus += mu(i) * mu(i); |
| 163 | + } |
| 164 | + double LowerBoundAndoLin = 1 + splus / sminus; |
| 165 | + UpdateBounds(LowerBoundAndoLin, "Ando/Lin lower bound"); |
| 166 | + std::cerr << "Ando Lin lower bound = " << LowerBoundAndoLin << "\n"; |
| 167 | + // |
| 168 | + double nPlus = 0; |
| 169 | + double nMinus = 0; |
| 170 | + for (int i = 0; i < n; i++) { |
| 171 | + if (mu(i) > 0) |
| 172 | + nPlus++; |
| 173 | + if (mu(i) < 0) |
| 174 | + nMinus++; |
| 175 | + } |
| 176 | + double LowerBoundInertia = 1 + std::max(nMinus / nPlus, nPlus / nMinus); |
| 177 | + UpdateBounds(LowerBoundInertia, "Inertia lower bound"); |
| 178 | + std::cerr << "Elphick Wocjan inertial lower bound = " << LowerBoundInertia |
| 179 | + << "\n"; |
| 180 | + // |
| 181 | + std::cerr << "Realizing method=" << RealizingBound |
| 182 | + << " BestLowerBound=" << BestLowerBound << "\n"; |
| 183 | + std::cerr << "Normal termination of ComputeChromaticLowerBounds\n"; |
| 184 | + } catch (TerminalException const &e) { |
| 185 | + exit(e.eVal); |
| 186 | + } |
| 187 | +} |
0 commit comments