1 #include "matrixconstructiondefault.h"
2 #include "../matrix/invert.h"
3 #include "../matrix/rank.h"
4 #include "../yal/logger.h"
5
6 using namespace sympol;
7 using namespace yal;
8
9 static LoggerPtr logger(Logger::getLogger("SymMatrixD"));
10
construct(const Polyhedron & poly)11 bool MatrixConstructionDefault::construct(const Polyhedron& poly) {
12 YALLOG_DEBUG(logger, "matrix construction with default");
13 const ulong matrixRank = poly.dimension() - 1;
14 const ulong matrixRows = poly.rows();
15 typedef matrix::Matrix<mpq_class> QMatrix;
16
17 // compute column defect of inequality matrix
18 QMatrix* Acopy = new QMatrix(poly.rows(), matrixRank);
19 {ulong j=0;
20 BOOST_FOREACH(const QArray& row, poly.rowPair()) {
21 for (ulong i=0; i<matrixRank; ++i) {
22 mpq_set(Acopy->at(j,i).get_mpq_t(), row[i+1]);
23 }
24 ++j;
25 }}
26 matrix::Rank<QMatrix> r(Acopy);
27 // freeColumns will contain the column indices we may skip to get a full dimensional polyhedron
28 std::set<uint> freeColumns;
29 r.columnDefect(std::inserter(freeColumns, freeColumns.end()));
30 delete Acopy;
31
32 typedef std::map<mpq_class,uint> WeightMap;
33 WeightMap weights;
34 mpq_t temp;
35 mpq_init(temp);
36
37 //
38 // calculate vertex weights
39 //
40 QMatrix Q(matrixRank - freeColumns.size());
41 ulong iQ = 0, jQ = 0;
42 for (ulong i=0; i<matrixRank; ++i) {
43 if (freeColumns.count(i)) {
44 YALLOG_DEBUG3(logger, "free column #" << (i+1));
45 continue;
46 }
47 jQ = 0;
48 for (ulong j=0; j<matrixRank; ++j) {
49 if (freeColumns.count(j))
50 continue;
51 BOOST_ASSERT(iQ < matrixRank - freeColumns.size());
52 BOOST_ASSERT(jQ < matrixRank - freeColumns.size());
53 BOOST_FOREACH(const QArray& row, poly.rowPair()) {
54 mpq_mul(temp, row[i+1], row[j+1]);
55 mpq_add(Q.at(iQ,jQ).get_mpq_t(), Q.at(iQ,jQ).get_mpq_t(), temp);
56 //Q.at(i,j) += row[i] * row[j];
57 }
58 ++jQ;
59 }
60 ++iQ;
61 }
62
63 YALLOG_DEBUG2(logger, "supposed rank of Q = " << matrixRank << "; rank defect = " << freeColumns.size() << "; Q = " << std::endl << Q);
64
65 QMatrix Qinv(Q.rows());
66 if (!matrix::Invert<QMatrix>(&Q).invert(&Qinv)) {
67 YALLOG_ERROR(logger, "could not invert matrix");
68 return false;
69 }
70
71 YALLOG_INFO(logger, "matrix inversion complete");
72 YALLOG_DEBUG3(logger, "Qinv = " << std::endl << Qinv);
73
74 m_zMatrix = new matrix::ZMatrix(matrixRows);
75
76 uint weightIndex = 0;
77 ulong i = 0, j = 0;
78 BOOST_FOREACH(const QArray& row1, poly.rowPair()) {
79 j = 0;
80 BOOST_FOREACH(const QArray& row2, poly.rowPair()) {
81 if (i < j)
82 break;
83 mpq_class newWeight;
84 ulong kQ = 0, lQ = 0;
85 for (uint k = 0; k < matrixRank; ++k) {
86 if (freeColumns.count(k))
87 continue;
88 lQ = 0;
89 for (uint l = 0; l < matrixRank; ++l) {
90 if (freeColumns.count(l))
91 continue;
92 BOOST_ASSERT(kQ < matrixRank - freeColumns.size());
93 BOOST_ASSERT(lQ < matrixRank - freeColumns.size());
94 mpq_mul(temp, row1[k+1], row2[l+1]);
95 mpq_mul(temp, temp, Qinv.at(kQ,lQ).get_mpq_t());
96 mpq_add(newWeight.get_mpq_t(), newWeight.get_mpq_t(), temp);
97 //newWeight += Qinv.at(k,l) * row1[k] * row2[l];
98 ++lQ;
99 }
100 ++kQ;
101 }
102 std::pair<WeightMap::iterator, bool> suc = weights.insert(std::make_pair(newWeight, weightIndex));
103 m_zMatrix->at(i,j) = suc.first->second;
104 m_zMatrix->at(j,i) = suc.first->second;
105 if (suc.second)
106 ++weightIndex;
107 ++j;
108 }
109 ++i;
110 }
111 mpq_clear(temp);
112
113 m_zMatrix->k() = weightIndex;
114
115 YALLOG_DEBUG(logger, "zMatrix = " << std::endl << *m_zMatrix);
116
117 initData(poly, weightIndex);
118 return true;
119 }
120
121
weightAt(unsigned int i,unsigned int j) const122 unsigned int MatrixConstructionDefault::weightAt(unsigned int i, unsigned int j) const {
123 BOOST_ASSERT(m_zMatrix != 0);
124 return m_zMatrix->at(i,j);
125 }
126