1 /* Frobby: Software for monomial ideal computations.
2    Copyright (C) 2009 Bjarke Hammersholt Roune (www.broune.com)
3 
4    This program is free software; you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation; either version 2 of the License, or
7    (at your option) any later version.
8 
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13 
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see http://www.gnu.org/licenses/.
16 */
17 #include "stdinc.h"
18 #include "SatBinomIdeal.h"
19 
20 #include "BigIdeal.h"
21 #include "Matrix.h"
22 
23 #include <sstream>
24 #include <limits>
25 
SatBinomIdeal()26 SatBinomIdeal::SatBinomIdeal() {
27 }
28 
SatBinomIdeal(const VarNames & names)29 SatBinomIdeal::SatBinomIdeal(const VarNames& names):
30   _names(names) {
31 }
32 
insert(const vector<mpz_class> & binom)33 void SatBinomIdeal::insert(const vector<mpz_class>& binom) {
34   ASSERT(binom.size() == getVarCount());
35 
36   _gens.push_back(binom);
37 }
38 
getGenerator(size_t index) const39 const vector<mpz_class>& SatBinomIdeal::getGenerator(size_t index) const {
40   ASSERT(index < getGeneratorCount());
41 
42   return _gens[index];
43 }
44 
getGeneratorCount() const45 size_t SatBinomIdeal::getGeneratorCount() const {
46   return _gens.size();
47 }
48 
print(FILE * out) const49 void SatBinomIdeal::print(FILE* out) const {
50   ostringstream tmp;
51   print(tmp);
52   fputs(tmp.str().c_str(), out);
53 }
54 
print(ostream & out) const55 void SatBinomIdeal::print(ostream& out) const {
56   out << "/---- SatBinomIdeal of " << _gens.size() << " generators:\n";
57   for (vector<vector<mpz_class> >::const_iterator it = _gens.begin();
58        it != _gens.end(); ++it) {
59     for (vector<mpz_class>::const_iterator entry = it->begin();
60          entry != it->end(); ++entry)
61       out << *entry << ' ';
62     out << '\n';
63   }
64   out << "----/ End of list.\n";
65 }
66 
clearAndSetNames(const VarNames & names)67 void SatBinomIdeal::clearAndSetNames(const VarNames& names) {
68   _gens.clear();
69   _names = names;
70 }
71 
getVarCount() const72 size_t SatBinomIdeal::getVarCount() const {
73   return _names.getVarCount();
74 }
75 
renameVars(const VarNames & names)76 void SatBinomIdeal::renameVars(const VarNames& names) {
77   ASSERT(names.getVarCount() == getVarCount());
78 
79   _names = names;
80 }
81 
getNames() const82 const VarNames& SatBinomIdeal::getNames() const {
83   return _names;
84 }
85 
clear()86 void SatBinomIdeal::clear() {
87   _gens.clear();
88   _names.clear();
89 }
90 
reserve(size_t size)91 void SatBinomIdeal::reserve(size_t size) {
92   _gens.reserve(size);
93 }
94 
removeGeneratorsWithLeadingZero()95 void SatBinomIdeal::removeGeneratorsWithLeadingZero() {
96   size_t gen = 0;
97   while (gen < getGeneratorCount()) {
98     if (getGenerator(gen)[0] == 0) {
99       _gens[gen] = _gens.back();
100       _gens.pop_back();
101     } else
102       ++gen;
103   }
104 }
105 
removeGeneratorsWithoutLeadingZero()106 void SatBinomIdeal::removeGeneratorsWithoutLeadingZero() {
107   size_t gen = 0;
108   while (gen < getGeneratorCount()) {
109     if (getGenerator(gen)[0] != 0) {
110       _gens[gen] = _gens.back();
111       _gens.pop_back();
112     } else
113       ++gen;
114   }
115 }
116 
getInitialIdeal(BigIdeal & ideal) const117 void SatBinomIdeal::getInitialIdeal(BigIdeal& ideal) const {
118   ideal.clearAndSetNames(getNames());
119   ideal.reserve(getGeneratorCount());
120 
121   for (size_t gen = 0; gen < getGeneratorCount(); ++gen) {
122     ideal.newLastTerm();
123     for (size_t var = 0; var < getVarCount(); ++var)
124       if (getGenerator(gen)[var] > 0)
125         ideal.getLastTermExponentRef(var) = getGenerator(gen)[var];
126   }
127 }
128 
getLastBinomRef()129 vector<mpz_class>& SatBinomIdeal::getLastBinomRef() {
130   ASSERT(!_gens.empty());
131 
132   return _gens.back();
133 }
134 
newLastTerm()135 void SatBinomIdeal::newLastTerm() {
136   _gens.resize(_gens.size() + 1);
137   _gens.back().resize(_names.getVarCount());
138 }
139 
hasZeroEntry() const140 bool SatBinomIdeal::hasZeroEntry() const {
141   for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
142     for (size_t var = 0; var < getVarCount(); ++var)
143       if (getGenerator(gen)[var] == 0)
144         return true;
145   return false;
146 }
147 
initialIdealIsWeaklyGeneric() const148 bool SatBinomIdeal::initialIdealIsWeaklyGeneric() const {
149   vector<mpz_class> v(getVarCount());
150 
151   for (size_t gen1 = 0; gen1 < getGeneratorCount(); ++gen1) {
152     for (size_t gen2 = gen1 + 1; gen2 < getGeneratorCount(); ++gen2) {
153       const vector<mpz_class>& g1 = getGenerator(gen1);
154       const vector<mpz_class>& g2 = getGenerator(gen2);
155 
156       // Skip if g1 and g2 are different in each entry.
157       bool sharesEntry = false;
158       for (size_t var = 0; var < getVarCount(); ++var) {
159         if (g1[var] == g2[var] && g1[var] > 0) {
160           sharesEntry = true;
161           break;
162         }
163       }
164       if (!sharesEntry)
165         continue;
166 
167       if (isPointFreeBody(g1, g2))
168         return false;
169     }
170   }
171 
172   return true;
173 }
174 
isPointFreeBody(const vector<mpz_class> & a,const vector<mpz_class> & b) const175 bool SatBinomIdeal::isPointFreeBody(const vector<mpz_class>& a,
176                                     const vector<mpz_class>& b) const {
177   ASSERT(a.size() == getVarCount());
178   ASSERT(b.size() == getVarCount());
179 
180   vector<mpz_class> rhs(getVarCount());
181 
182   // Set rhs to max(0,g1,g2)-1.
183   for (size_t var = 0; var < getVarCount(); ++var) {
184     rhs[var] = a[var] > b[var] ? a[var] : b[var];
185     if (rhs[var] < 0)
186       rhs[var] = 0;
187     rhs[var] -= 1;
188   }
189 
190   return !isDominating(rhs);
191 }
192 
isPointFreeBody(const vector<mpz_class> & a,const vector<mpz_class> & b,const vector<mpz_class> & c) const193 bool SatBinomIdeal::isPointFreeBody(const vector<mpz_class>& a,
194                                     const vector<mpz_class>& b,
195 									const vector<mpz_class>& c) const {
196   ASSERT(a.size() == getVarCount());
197   ASSERT(b.size() == getVarCount());
198 
199   vector<mpz_class> rhs(getVarCount());
200 
201   // Set rhs to max(0,g1,g1+g2)-1.
202   for (size_t var = 0; var < getVarCount(); ++var) {
203     rhs[var] = a[var] > b[var] ? a[var] : b[var];
204 	rhs[var] = rhs[var] > c[var] ? rhs[var] : c[var];
205     if (rhs[var] < 0)
206       rhs[var] = 0;
207     rhs[var] -= 1;
208   }
209 
210   return !isDominating(rhs);
211 }
212 
isInterior(const vector<mpz_class> & a,const vector<mpz_class> & b) const213 bool SatBinomIdeal::isInterior(const vector<mpz_class>& a,
214 							   const vector<mpz_class>& b) const {
215   ASSERT(a.size() == b.size());
216 
217   if (!isPointFreeBody(a, b))
218 	return false;
219 
220   for (size_t var = 1; var < a.size(); ++var)
221 	if (a[var] <= 0 && b[var] <= 0)
222 	  return false;
223   return true;
224 }
225 
226 namespace {
hasCycle(size_t gen,vector<char> & color,const SatBinomIdeal & ideal)227   bool hasCycle(size_t gen, vector<char>& color, const SatBinomIdeal& ideal) {
228 	// 0 = Not seen before.
229 	// 1 = Exploring now.
230 	// 2 = Already explored. Does not lead to cycle.
231 	if (color[gen] == 1)
232 	  return true;
233 	if (color[gen] == 2)
234 	  return false;
235 	color[gen] = 1;
236 	for (size_t g = 0; g < ideal.getGeneratorCount(); ++g)
237 	  if (ideal.isInteriorEdge(gen, g) &&
238 		  !ideal.isTerminatingEdge(gen, g) &&
239 		  hasCycle(g, color, ideal))
240 		return true;
241 	color[gen] = 2;
242 	return false;
243   }
244 }
245 
validate() const246 bool SatBinomIdeal::validate() const {
247   // check the graph satisifies what we think it should.
248 
249   bool generic = !hasZeroEntry();
250   if (!generic)
251 	return true;
252 
253   // termination.
254   vector<char> color(getGeneratorCount());
255   for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
256 	if (hasCycle(gen, color, *this))
257 	  return false;
258 
259   // out-degree 1 when generic
260   for (size_t from = 0; from < getGeneratorCount(); ++from) {
261 	const vector<mpz_class>& fromGen = getGenerator(from);
262 
263 	size_t outDegree = 0;
264 	for (size_t to = 0; to < getGeneratorCount(); ++to) {
265 	  if (isInteriorEdge(from, to))
266 		++outDegree;
267 	}
268 	if (isInterior(fromGen, fromGen)) {
269 	  if (outDegree != 0)
270 		return false;
271 	} else {
272 	  if (outDegree != 1)
273 		return false;
274 	}
275   }
276 
277   return true;
278 }
279 
isInteriorEdge(size_t from,size_t to) const280 bool SatBinomIdeal::isInteriorEdge(size_t from, size_t to) const {
281   const vector<mpz_class>& fromGen = getGenerator(from);
282   const vector<mpz_class>& toGen = getGenerator(to);
283 
284   if (isInterior(fromGen, fromGen))
285 	return false;
286   if (isInterior(toGen, toGen))
287 	return false;
288 
289   vector<mpz_class> sum(fromGen.size());
290   for (size_t var = 0; var < fromGen.size(); ++var)
291 	sum[var] = fromGen[var] + toGen[var];
292   return isInterior(toGen, sum);
293 }
294 
isTerminatingEdge(size_t from,size_t to) const295 bool SatBinomIdeal::isTerminatingEdge(size_t from, size_t to) const {
296   if (!isInteriorEdge(from, to))
297 	return false;
298 
299   const vector<mpz_class> fromGen = getGenerator(from);
300   const vector<mpz_class> toGen = getGenerator(to);
301 
302   vector<mpz_class> sum(fromGen.size());
303   for (size_t var = 0; var < fromGen.size(); ++var)
304 	sum[var] = fromGen[var] + toGen[var];
305   return isPointFreeBody(fromGen, sum);
306 }
307 
getDoubleTriangleCount(mpz_class & count) const308 void SatBinomIdeal::getDoubleTriangleCount(mpz_class& count) const {
309   vector<mpz_class> sum(getVarCount());
310 
311   count = 0;
312   for (size_t gen1 = 0; gen1 < getGeneratorCount(); ++gen1) {
313     for (size_t gen2 = gen1 + 1; gen2 < getGeneratorCount(); ++gen2) {
314       const vector<mpz_class>& g1 = getGenerator(gen1);
315       const vector<mpz_class>& g2 = getGenerator(gen2);
316 
317       // Set sum = g1 + g2.
318       for (size_t var = 0; var < getVarCount(); ++var)
319         sum[var] = g1[var] + g2[var];
320 
321       if (isPointFreeBody(g1, sum) && isPointFreeBody(g2, sum))
322         ++count;
323     }
324   }
325 }
326 
isGeneric() const327 bool SatBinomIdeal::isGeneric() const {
328   return !hasZeroEntry() && initialIdealIsWeaklyGeneric();
329 }
330 
operator =(const SatBinomIdeal & ideal)331 SatBinomIdeal& SatBinomIdeal::operator=(const SatBinomIdeal& ideal) {
332   _gens = ideal._gens;
333   _names = ideal._names;
334   return *this;
335 }
336 
isDominating(const vector<mpz_class> & v) const337 bool SatBinomIdeal::isDominating(const vector<mpz_class>& v) const {
338   for (size_t gen = 0; gen < getGeneratorCount(); ++gen) {
339     bool dom = true;
340     for (size_t var = 0; var < getVarCount(); ++var) {
341       if (v[var] < getGenerator(gen)[var]) {
342         dom = false;
343         break;
344       }
345     }
346     if (dom)
347       return true;
348   }
349   return false;
350 }
351 
isGenerator(const vector<mpz_class> & v) const352 bool SatBinomIdeal::isGenerator
353 (const vector<mpz_class>& v) const {
354   for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
355     if (getGenerator(gen) == v)
356       return true;
357   return false;
358 }
359 
projectVar(size_t var)360 void SatBinomIdeal::projectVar(size_t var) {
361   ASSERT(var < getVarCount());
362 
363   for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
364     _gens[gen].erase(_gens[gen].begin() + var);
365   _names.projectVar(var);
366 }
367 
getMatrix(Matrix & matrix) const368 void SatBinomIdeal::getMatrix(Matrix& matrix) const {
369   matrix.resize(getGeneratorCount(), getVarCount());
370   for (size_t gen = 0; gen < getGeneratorCount(); ++gen)
371 	for (size_t var = 0; var < getVarCount(); ++var)
372 	  matrix(gen, var) = _gens[gen][var];
373 }
374