Frobby  0.9.0
SliceFacade.cpp
Go to the documentation of this file.
1 /* Frobby: Software for monomial ideal computations.
2  Copyright (C) 2007 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 "SliceFacade.h"
19 
20 #include "BigTermConsumer.h"
21 #include "CoefBigTermConsumer.h"
22 #include "TermTranslator.h"
23 #include "BigIdeal.h"
24 #include "Ideal.h"
25 #include "Term.h"
26 #include "MsmStrategy.h"
29 #include "DebugStrategy.h"
30 #include "DecomRecorder.h"
31 #include "TermGrader.h"
32 #include "OptimizeStrategy.h"
34 #include "HilbertStrategy.h"
35 #include "IOHandler.h"
36 #include "BigPolynomial.h"
38 #include "CoefBigTermRecorder.h"
39 #include "CanonicalTermConsumer.h"
40 #include "VarSorter.h"
41 #include "StatisticsStrategy.h"
43 #include "SizeMaxIndepSetAlg.h"
44 #include "SliceParams.h"
45 #include "error.h"
46 #include "display.h"
47 
48 #include <iterator>
49 
50 SliceFacade::SliceFacade(const SliceParams& params, const DataType& output):
51  Facade(params.getPrintActions()),
52  _params(params) {
53  _split = SplitStrategy::createStrategy(params.getSplit().c_str());
54  _common.readIdealAndSetOutput(params, output);
55 }
56 
58  const BigIdeal& ideal,
59  BigTermConsumer& consumer):
60  Facade(params.getPrintActions()),
61  _params(params) {
62  _split = SplitStrategy::createStrategy(params.getSplit().c_str());
63  _common.setIdealAndIdealOutput(params, ideal, consumer);
64 }
65 
67  const BigIdeal& ideal,
68  CoefBigTermConsumer& consumer):
69  Facade(params.getPrintActions()),
70  _params(params) {
71  _split = SplitStrategy::createStrategy(params.getSplit().c_str());
72  _common.setIdealAndPolyOutput(params, ideal, consumer);
73 }
74 
76 }
77 
80  beginAction("Computing multigraded Hilbert-Poincare series.");
81 
82  auto_ptr<CoefTermConsumer> consumer = _common.makeTranslatedPolyConsumer();
83 
84  consumer->consumeRing(_common.getNames());
85  consumer->beginConsuming();
86  HilbertStrategy strategy(consumer.get(), _split.get());
88  consumer->doneConsuming();
89 
90  endAction();
91 }
92 
95  beginAction("Computing univariate Hilbert-Poincare series.");
96 
97  auto_ptr<CoefTermConsumer> consumer =
99 
100  consumer->consumeRing(_common.getNames());
101  consumer->beginConsuming();
102  HilbertStrategy strategy(consumer.get(), _split.get());
104  consumer->doneConsuming();
105 
106  endAction();
107 }
108 
112 }
113 
116 
118  if (codimension) {
119  // convert to mpz_class before increment to ensure no overflow.
120  return mpz_class(_common.getIdeal().getVarCount()) + 1;
121  } else
122  return -1;
123  }
124 
125  // todo: inline?
126  takeRadical();
127 
128  beginAction("Preparing to compute dimension.");
129 
130  vector<mpz_class> v;
131  fill_n(back_inserter(v), _common.getIdeal().getVarCount(), -1);
132 
133  endAction();
134 
135  mpz_class minusCodimension;
136 #ifdef DEBUG
137  // Only define hasComponents when DEBUG is defined since otherwise
138  // GCC will warn about hasComponents not being used.
139  bool hasComponents =
140 #endif
141  solveIrreducibleDecompositionProgram(v, minusCodimension, false);
142  ASSERT(hasComponents);
143 
144  if (codimension)
145  return -minusCodimension;
146  else
147  return v.size() + minusCodimension;
148 }
149 
152 
153  size_t varCount = _common.getIdeal().getVarCount();
154 
155  Ideal irreducibleDecom(varCount);
156  {
157  DecomRecorder recorder(&irreducibleDecom);
158  produceEncodedIrrDecom(recorder);
159  }
160 
162  ("Computing primary decomposition from irreducible decomposition.");
163 
164  // Do intersection of each component also using irreducible
165  // decomposition of the dual. We can't use the Alexander dual
166  // methods, since those switch around the translator to emit altered
167  // big integers, while keeping the small integers the same, but we
168  // want to keep this in small integers. So we have to do the dual
169  // thing here.
170 
171  // To get actual supports.
172  _common.getTranslator().setInfinityPowersToZero(irreducibleDecom);
173 
174  // To collect same-support vectors together.
175  irreducibleDecom.sortReverseLex();
176 
177  Term lcm(varCount);
178  irreducibleDecom.getLcm(lcm);
179 
180  Term tmp(varCount);
181  Term support(varCount);
182 
183  _common.getIdeal().clear();
184  Ideal& primaryComponentDual = _common.getIdeal();
185  Ideal primaryComponent(varCount);
186 
187  DecomRecorder recorder(&primaryComponent);
188 
189  auto_ptr<TermConsumer> consumer = _common.makeTranslatedIdealConsumer();
190  consumer->consumeRing(_common.getNames());
191  consumer->beginConsumingList();
192 
193  Ideal::const_iterator stop = irreducibleDecom.end();
194  Ideal::const_iterator it = irreducibleDecom.begin();
195  while (it != stop) {
196  // Get all vectors with same support.
197  support = *it;
198  do {
199  tmp.encodedDual(*it, lcm);
200  primaryComponentDual.insert(tmp);
201  ++it;
202  } while (it != stop && support.hasSameSupport(*it));
203  ASSERT(!primaryComponentDual.isZeroIdeal());
204 
205  _common.getTranslator().addPurePowersAtInfinity(primaryComponentDual);
206  {
207  MsmStrategy strategy(&recorder, _split.get());
209  }
210  _common.getTranslator().setInfinityPowersToZero(primaryComponent);
211 
212  consumer->beginConsuming();
213  for (Ideal::const_iterator dualTerm = primaryComponent.begin();
214  dualTerm != primaryComponent.end(); ++dualTerm) {
215  tmp.encodedDual(*dualTerm, lcm);
216  consumer->consume(tmp);
217  }
218  consumer->doneConsuming();
219 
220  primaryComponent.clear();
221  primaryComponentDual.clear();
222  }
223 
224  consumer->doneConsumingList();
225 
226  endAction();
227 }
228 
231  beginAction("Computing maximal staircase monomials.");
232 
233  auto_ptr<TermConsumer> consumer = _common.makeTranslatedIdealConsumer();
234  consumer->consumeRing(_common.getNames());
235  MsmStrategy strategy(consumer.get(), _split.get());
237 
238  endAction();
239 }
240 
243 
244  beginAction("Preparing to compute maximal standard monomials.");
246  endAction();
248 }
249 
250 void SliceFacade::computeAlexanderDual(const vector<mpz_class>& point) {
252  ASSERT(point.size() == _common.getIdeal().getVarCount());
253 
254  beginAction("Ensuring specified point is divisible by lcm.");
255  vector<mpz_class> lcm;
256  getLcmOfIdeal(lcm);
257 
258  for (size_t var = 0; var < lcm.size(); ++var) {
259  if (lcm[var] > point[var]) {
260  endAction();
262  ("The specified point to dualize on is not divisible by the "
263  "least common multiple of the minimal generators of the ideal.");
264  }
265  }
266  endAction();
267 
268  beginAction("Preparing to compute Alexander dual.");
269  _common.getTranslator().dualize(point);
270  endAction();
271 
273 }
274 
277 
278  beginAction("Computing lcm for Alexander dual.");
279  vector<mpz_class> lcm;
280  getLcmOfIdeal(lcm);
281  endAction();
282 
284 }
285 
288 
289  size_t varCount = _common.getIdeal().getVarCount();
290 
291  // Obtain generators of radical from irreducible decomposition.
292  Ideal radical(varCount);
293  {
294  Ideal decom(varCount);
295  DecomRecorder recorder(&decom);
296  produceEncodedIrrDecom(recorder);
297 
298  beginAction("Computing associated primes from irreducible decomposition.");
299 
300  Term tmp(varCount);
301  Ideal::const_iterator stop = decom.end();
302  for (Ideal::const_iterator it = decom.begin(); it != stop; ++it) {
303  for (size_t var = 0; var < varCount; ++var) {
304  // We cannot just check whether (*it)[var] == 0, since the
305  // added fake pure powers map to zero but are not themselves
306  // zero.
307  if (_common.getTranslator().getExponent(var, (*it)[var]) == 0)
308  tmp[var] = 0;
309  else
310  tmp[var] = 1;
311  }
312  radical.insert(tmp);
313  }
314  }
315 
316  radical.removeDuplicates();
317 
318 
319  // Output associated primes.
321  auto_ptr<TermConsumer> consumer = _common.makeTranslatedIdealConsumer();
322 
323  consumer->consumeRing(_common.getNames());
324  consumer->beginConsuming();
325  Term tmp(varCount);
326  Ideal::const_iterator stop = radical.end();
327  for (Ideal::const_iterator it = radical.begin(); it != stop; ++it) {
328  tmp = *it;
329  consumer->consume(tmp);
330  }
331  consumer->doneConsuming();
332 
333  endAction();
334 }
335 
337 (const vector<mpz_class>& grading,
338  mpz_class& optimalValue,
339  bool reportAllSolutions) {
341  ASSERT(grading.size() == _common.getIdeal().getVarCount());
342 
343  beginAction("Preparing to solve optimization program.");
346  endAction();
347 
348  return solveProgram(grading, optimalValue, reportAllSolutions);
349 }
350 
352 (const vector<mpz_class>& grading,
353  mpz_class& optimalValue,
354  bool reportAllSolutions) {
356  ASSERT(grading.size() == _common.getIdeal().getVarCount());
357 
359  return solveProgram(grading, optimalValue, reportAllSolutions);
360 }
361 
364  beginAction("Computing irreducible decomposition.");
365 
367  MsmStrategy strategy(&consumer, _split.get());
368 
369  consumer.consumeRing(_common.getNames());
371 
372  endAction();
373 }
374 
375 bool SliceFacade::solveProgram(const vector<mpz_class>& grading,
376  mpz_class& optimalValue,
377  bool reportAllSolutions) {
379  ASSERT(grading.size() == _common.getIdeal().getVarCount());
380 
383  ("Turning off Independence splits as they are not supported\n"
384  "for optimization.");
386  }
387 
391  ("Bound simplification requires using the bound to eliminate\n"
392  "non-improving slices, which has been turned off. Am now turning\n"
393  "this on.");
395  }
396 
397  beginAction("Solving optimization program.");
398 
399  OptimizeStrategy::BoundSetting boundSetting;
403  } else if (_params.getUseBoundElimination())
405  else
406  boundSetting = OptimizeStrategy::DoNotUseBound;
407 
408  TermGrader grader(grading, _common.getTranslator());
409  OptimizeStrategy strategy
410  (grader, _split.get(), reportAllSolutions, boundSetting);
412 
413  endAction();
414 
415  const Ideal& solution = strategy.getMaximalSolutions();
416 
417  auto_ptr<TermConsumer> consumer = _common.makeTranslatedIdealConsumer();
418  consumer->consumeRing(_common.getNames());
419  consumer->consume(solution);
420 
421  if (solution.isZeroIdeal())
422  return false;
423  else {
424  optimalValue = strategy.getMaximalValue();
425  return true;
426  }
427 }
428 
430  return _common.hasIdeal();
431 }
432 
435 
436  beginAction("Taking radical of ideal.");
437 
438  bool skip = false;
441  if (lcm.isSquareFree())
442  skip = true;
443 
444  if (!skip) {
448  }
449 
451 
452  endAction();
453 }
454 
455 void SliceFacade::getLcmOfIdeal(vector<mpz_class>& bigLcm) {
457 
460 
461  bigLcm.clear();
462  bigLcm.reserve(_common.getIdeal().getVarCount());
463  for (size_t var = 0; var < _common.getIdeal().getVarCount(); ++var)
464  bigLcm.push_back(_common.getTranslator().getExponent(var, lcm));
465 }
466 
471 
472  SliceStrategy* strategyWithOptions = &strategy;
473 
474  auto_ptr<SliceStrategy> debugStrategy;
475  if (_params.getPrintDebug()) {
476  debugStrategy.reset
477  (new DebugStrategy(strategyWithOptions, stderr));
478  strategyWithOptions = debugStrategy.get();
479  }
480 
481  auto_ptr<SliceStrategy> statisticsStrategy;
482  if (_params.getPrintStatistics()) {
483  statisticsStrategy.reset
484  (new StatisticsStrategy(strategyWithOptions, stderr));
485  strategyWithOptions = statisticsStrategy.get();
486  }
487 
488  ASSERT(strategyWithOptions != 0);
489  strategyWithOptions->run(_common.getIdeal());
490 }
void setIdealAndPolyOutput(const CommonParams &params, const BigIdeal &input, CoefBigTermConsumer &output)
Use given ideal and support polynomial output.
Cont::const_iterator const_iterator
Definition: Ideal.h:43
Eliminate non-improving slices, achieving a branch-and-bound algorithm in place of the usual backtrac...
auto_ptr< CoefTermConsumer > makeToUnivariatePolyConsumer()
static void encodedDual(Exponent *res, const Exponent *dualOf, const Exponent *point, size_t varCount)
The parameter dualOf is interpreted to encode an irreducible ideal, and the dual of that reflected in...
Definition: Term.h:497
auto_ptr< CoefTermConsumer > makeTranslatedPolyConsumer()
bool containsIdentity() const
Definition: Ideal.cpp:65
void computeAlexanderDual()
Compute the Alexander dual of the ideal.
This is the super class of all facades.
Definition: Facade.h:32
void setInfinityPowersToZero(Ideal &ideal) const
The method addPurePowersAtInfinity adds high exponents that map to zero.
bool getUseSimplification() const
Apply simplification to the state of the algorithm when possible.
virtual void run(const Ideal &ideal)=0
Run the Slice algorithm.
auto_ptr< SplitStrategy > _split
Definition: SliceFacade.h:225
void takeRadical()
OptimizeStrategy optimizes a function on the maximal standard monomials of a monomial ideal using bra...
#define ASSERT(X)
Definition: stdinc.h:85
void minimize()
Definition: Ideal.cpp:501
void getLcmOfIdeal(vector< mpz_class > &lcm)
const_iterator begin() const
Definition: Ideal.h:48
void computeMultigradedHilbertSeries()
Compute the numerator of the multigraded Hilbert-Poincare series.
Definition: SliceFacade.cpp:78
void sortReverseLex()
Definition: Ideal.cpp:510
Represents a monomial ideal with int exponents.
Definition: Ideal.h:27
const_iterator end() const
Definition: Ideal.h:49
void lcm(Word *res, const Word *resEnd, const Word *a, const Word *b)
bool isFirstComputation() const
A wrapper for a SliceStrategy that collects statistics on what is going on, while delegating everythi...
The intention of this class is to describe the different kinds of mathematical structures that Frobby...
Definition: DataType.h:29
bool getPrintStatistics() const
Returns whether to print statistics on what the algorithm did to standard error after it has run...
Definition: CommonParams.h:66
bool isZeroIdeal() const
Definition: Ideal.cpp:86
auto_ptr< TermConsumer > makeTranslatedIdealConsumer(bool split=false)
void codimension(const Ideal &ideal, mpz_t codim)
Compute the codimension of a monomial ideal.
Definition: frobby.cpp:437
A TermGrader assigns a value, the degree, to each monomial.
Definition: TermGrader.h:27
This class is used to transfer terms one at a time from one part of the program to another...
Definition: TermConsumer.h:36
const mpz_class & getExponent(size_t variable, Exponent exponent) const
This method translates from IDs to arbitrary precision integers.
void computeMaximalStandardMonomials()
Compute the maximal standard monomials of the ideal.
This file contains functions for printing strings to standard error.
virtual void consumeRing(const VarNames &names)
Tell the consumer which ring is being used.
void produceEncodedIrrDecom(TermConsumer &consumer)
const string & getSplit() const
Definition: SliceParams.h:30
void setIdealAndIdealOutput(const CommonParams &params, const BigIdeal &input, BigTermConsumer &output)
Use given ideal and support ideal output.
bool solveStandardProgram(const vector< mpz_class > &grading, mpz_class &value, bool reportAllSolutions)
Solve an optimization program over maximal standard monomials.
void endAction()
Prints to standard error the time since the last call to beginAction.
Definition: Facade.cpp:51
void displayNote(const string &msg)
Display msg to standard error in a way that indicates that this is something that the user should tak...
Definition: display.cpp:135
virtual void setUseIndependence(bool use)=0
This method should only be called before calling run().
bool getUseBoundSimplification() const
Returns whether to simplify slices by seeking to generate non-improving slices that are then eliminat...
Definition: SliceParams.h:44
static auto_ptr< SplitStrategy > createStrategy(const string &prefix)
Returns the strategy whose name has the given prefix.
void decrement()
Replaces var^v by var^(v-1).
const VarNames & getNames()
void readIdealAndSetOutput(const CommonParams &params, const DataType &output)
Read input ideal and support specified kind of output.
virtual void setUseSimplification(bool use)=0
This method should only be called before calling run().
SliceFacade(const SliceParams &params, const DataType &output)
Definition: SliceFacade.cpp:50
void computeIrreducibleDecomposition(bool encode)
Compute the unique irredundant set of irreducible ideals whose intersection equals ideal...
Make no use of the bound.
void removeDuplicates()
Definition: Ideal.cpp:634
size_t getVarCount() const
Definition: Ideal.h:56
void beginAction(const char *message)
Prints message to standard error if printing is turned on, and records the time when the action start...
Definition: Facade.cpp:38
mpz_class computeDimension(bool codimension=false)
Compute the Krull dimension of ideal.
void computeAssociatedPrimes()
Compute the associated primes of the ideal.
void useIndependenceSplits(bool value)
Definition: SliceParams.h:34
void dualize(const vector< mpz_class > &a)
Replaces var^v by var^(a[i] - v) except that var^0 is left alone.
TermTranslator & getTranslator()
SliceParams _params
Definition: SliceFacade.h:223
void useBoundElimination(bool value)
Definition: SliceParams.h:39
void setToZeroOne(TermTranslator &translator)
void addPurePowersAtInfinity(Ideal &ideal) const
Adds a generator of the form v^e, e > 0, for any variable v where generator of that form is not alrea...
This class describes the interface of a strategy object for the Slice Algorithm.
Definition: SliceStrategy.h:33
CommonParamsHelper _common
Definition: SliceFacade.h:224
Eliminate non-improving slices and simplify slices by trying to generate non-improving slices that ar...
BoundSetting
The values of BoundSetting indicate how to use the bound.
void clear()
Definition: Ideal.cpp:641
bool getUseBoundElimination() const
Returns whether to use branch-and-bound to speed up Slice optimization computations by eliminating no...
Definition: SliceParams.h:38
void computeMaximalStaircaseMonomials()
Compute the maximal staircase monomials of the ideal.
bool solveIrreducibleDecompositionProgram(const vector< mpz_class > &grading, mpz_class &optimalValue, bool reportAllSolutions)
Solve an optimization program over irreducible components.
void reportError(const string &errorMsg)
Definition: error.cpp:23
void insert(const Exponent *term)
Definition: Ideal.cpp:455
bool solveProgram(const vector< mpz_class > &grading, mpz_class &optimalValue, bool reportAllSolutions)
void computeUnivariateHilbertSeries()
Compute the numerator of the univariate Hilbert-Poincare series.
Definition: SliceFacade.cpp:93
void takeRadicalNoMinimize()
Replaces all generators with their support and does not remove any non-minimal generators this may pr...
Definition: Ideal.cpp:660
static bool hasSameSupport(const Exponent *a, const Exponent *b, size_t varCount)
Returns whether for every variable .
Definition: Term.h:438
void runSliceAlgorithmWithOptions(SliceStrategy &strategy)
void computePrimaryDecomposition()
Compute the unique "nicest" primary decomposition of the ideal.
Term represents a product of variables which does not include a coefficient.
Definition: Term.h:49
bool getUseIndependenceSplits() const
Definition: SliceParams.h:33
bool getPrintDebug() const
Returns whether to print information about what the algorithm is doing to standard error as it runs...
Definition: CommonParams.h:61
void getLcm(Exponent *lcm) const
Sets lcm to the least common multiple of all generators.
Definition: Ideal.cpp:157