1 | /* |
---|
2 | gfan.h Interface to gfan.cc |
---|
3 | |
---|
4 | $Author: monerjan $ |
---|
5 | $Date: 2009/11/03 06:57:32 $ |
---|
6 | $Header: /usr/local/Singular/cvsroot/kernel/gfan.h,v 1.13 2009/11/03 06:57:32 monerjan Exp $ |
---|
7 | $Id$ |
---|
8 | */ |
---|
9 | #ifdef HAVE_GFAN |
---|
10 | |
---|
11 | #ifndef GFAN_H |
---|
12 | #define GFAN_H |
---|
13 | |
---|
14 | #include "intvec.h" |
---|
15 | #include "intvec.h" |
---|
16 | |
---|
17 | #define p800 |
---|
18 | #ifdef p800 |
---|
19 | #include "../../cddlib/include/setoper.h" |
---|
20 | #include "../../cddlib/include/cdd.h" |
---|
21 | #include "../../cddlib/include/cddmp.h" |
---|
22 | #endif |
---|
23 | extern int gfanHeuristic; |
---|
24 | // extern dd_MatrixPtr ddLinealitySpace; |
---|
25 | |
---|
26 | // #ifdef gfanp |
---|
27 | // extern static float time_getConeNormals; |
---|
28 | // extern static float time_getCodim2Normals; |
---|
29 | // extern static float time_flip; |
---|
30 | // extern static float time_enqueue; |
---|
31 | // extern static float time_computeInv; |
---|
32 | // #endif |
---|
33 | //ideal getGB(ideal inputIdeal); |
---|
34 | // ideal gfan(ideal inputIdeal, int heuristic); |
---|
35 | lists gfan(ideal inputIdeal, int heuristic); |
---|
36 | //int dotProduct(intvec a, intvec b); |
---|
37 | //bool isParallel(intvec a, intvec b); |
---|
38 | |
---|
39 | class facet |
---|
40 | { |
---|
41 | private: |
---|
42 | /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */ |
---|
43 | intvec *fNormal; |
---|
44 | |
---|
45 | /** \brief An interior point of the facet*/ |
---|
46 | intvec *interiorPoint; |
---|
47 | |
---|
48 | /** \brief Universal Cone Number |
---|
49 | * The number of the cone the facet belongs to, Set in getConeNormals() |
---|
50 | */ |
---|
51 | int UCN; |
---|
52 | |
---|
53 | /** \brief The codim of the facet |
---|
54 | */ |
---|
55 | int codim; |
---|
56 | |
---|
57 | /** \brief The Groebner basis on the other side of a shared facet |
---|
58 | * |
---|
59 | * In order not to have to compute the flipped GB twice we store the basis we already get |
---|
60 | * when identifying search facets. Thus in the next step of the reverse search we can |
---|
61 | * just copy the old cone and update the facet and the gcBasis. |
---|
62 | * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB |
---|
63 | */ |
---|
64 | ideal flipGB; //The Groebner Basis on the other side, computed via gcone::flip |
---|
65 | |
---|
66 | public: |
---|
67 | /** \brief Boolean value to indicate whether a facet is flippable or not |
---|
68 | * This is also used to mark facets that nominally are flippable but which do |
---|
69 | * not intersect with the positive orthant. This check is done in gcone::getCodim2Normals |
---|
70 | */ |
---|
71 | bool isFlippable; //**flippable facet? */ |
---|
72 | bool isIncoming; //Is the facet incoming or outgoing in the reverse search? |
---|
73 | facet *next; //Pointer to next facet |
---|
74 | facet *prev; //Pointer to predecessor. Needed for the SearchList in noRevS |
---|
75 | facet *codim2Ptr; //Pointer to (codim-2)-facet. Bit of recursion here ;-) |
---|
76 | int numCodim2Facets; //#of (codim-2)-facets of this facet. Set in getCodim2Normals() |
---|
77 | unsigned numRays; //Number of spanning rays of the cone |
---|
78 | ring flipRing; //the ring on the other side of the facet |
---|
79 | |
---|
80 | /** The default constructor. */ |
---|
81 | facet(); |
---|
82 | /** Constructor for lower dimensional faces*/ |
---|
83 | facet(const int &n); |
---|
84 | /** The copy constructor */ |
---|
85 | facet(const facet& f); |
---|
86 | facet* shallowCopy(const facet& f); |
---|
87 | void shallowDelete(); |
---|
88 | /** The default destructor */ |
---|
89 | ~facet(); |
---|
90 | /** Comparison operator*/ |
---|
91 | // inline bool operator==(const facet *f,const facet *g); |
---|
92 | /** \brief Comparison of facets*/ |
---|
93 | inline bool areEqual(facet *f, facet *g); |
---|
94 | /** Stores the facet normal \param intvec*/ |
---|
95 | inline void setFacetNormal(intvec *iv); |
---|
96 | /** Hopefully returns the facet normal */ |
---|
97 | inline intvec *getFacetNormal(); |
---|
98 | /** Return a reference to the facet normal*/ |
---|
99 | inline const intvec *getRef2FacetNormal(); |
---|
100 | /** Method to print the facet normal*/ |
---|
101 | inline void printNormal(); |
---|
102 | /** Store the flipped GB*/ |
---|
103 | inline void setFlipGB(ideal I); |
---|
104 | /** Return the flipped GB*/ |
---|
105 | inline ideal getFlipGB(); |
---|
106 | /** Print the flipped GB*/ |
---|
107 | inline void printFlipGB(); |
---|
108 | /** Set the UCN */ |
---|
109 | inline void setUCN(int n); |
---|
110 | /** \brief Get the UCN |
---|
111 | * Returns the UCN iff this != NULL, else -1 |
---|
112 | */ |
---|
113 | inline int getUCN(); |
---|
114 | /** Store an interior point of the facet */ |
---|
115 | inline void setInteriorPoint(intvec *iv); |
---|
116 | inline intvec *getInteriorPoint(); |
---|
117 | inline const intvec *getRef2InteriorPoint(); |
---|
118 | /** \brief Debugging function |
---|
119 | * prints the facet normal an all (codim-2)-facets that belong to it |
---|
120 | */ |
---|
121 | volatile void fDebugPrint(); |
---|
122 | friend class gcone; |
---|
123 | }; |
---|
124 | |
---|
125 | /** |
---|
126 | *\brief Implements the cone structure |
---|
127 | * |
---|
128 | * A cone is represented by a linked list of facet normals |
---|
129 | * @see facet |
---|
130 | */ |
---|
131 | |
---|
132 | class gcone |
---|
133 | { |
---|
134 | private: |
---|
135 | // ring rootRing; //good to know this -> generic walk |
---|
136 | ideal inputIdeal; //the original |
---|
137 | ring baseRing; //the basering of the cone |
---|
138 | intvec *ivIntPt; //an interior point of the cone |
---|
139 | int UCN; //unique number of the cone |
---|
140 | int pred; //UCN of the cone this one is derived from |
---|
141 | static int counter; |
---|
142 | |
---|
143 | public: |
---|
144 | /** \brief Pointer to the first facet */ |
---|
145 | facet *facetPtr; //Will hold the adress of the first facet; set by gcone::getConeNormals |
---|
146 | #ifdef gfanp |
---|
147 | static float time_getConeNormals; |
---|
148 | static float time_getCodim2Normals; |
---|
149 | static float t_getExtremalRays; |
---|
150 | static float t_ddPolyh; |
---|
151 | static float time_flip; |
---|
152 | static float time_flip2; |
---|
153 | static float t_areEqual; |
---|
154 | static float t_ffG; |
---|
155 | static float t_markings; |
---|
156 | static float t_dd; |
---|
157 | static float t_kStd; |
---|
158 | static float time_enqueue; |
---|
159 | static float time_computeInv; |
---|
160 | static float t_ddMC; |
---|
161 | static float t_mI; |
---|
162 | static float t_iP; |
---|
163 | static float t_isParallel; |
---|
164 | static unsigned parallelButNotEqual; |
---|
165 | static unsigned numberOfFacetChecks; |
---|
166 | #endif |
---|
167 | /** Matrix to contain the homogeneity/lineality space */ |
---|
168 | static dd_MatrixPtr dd_LinealitySpace; |
---|
169 | static int lengthOfSearchList; |
---|
170 | /** Maximum size of the searchlist*/ |
---|
171 | static int maxSize; |
---|
172 | /** is the ideal homogeneous? */ |
---|
173 | static bool hasHomInput; |
---|
174 | /** # of variables in the ring */ |
---|
175 | int numVars; //#of variables in the ring |
---|
176 | /** The hilbert function - for the homogeneous case*/ |
---|
177 | static intvec *hilbertFunction; |
---|
178 | |
---|
179 | /** # of facets of the cone |
---|
180 | * This value is set by gcone::getConeNormals |
---|
181 | */ |
---|
182 | int numFacets; //#of facets of the cone |
---|
183 | |
---|
184 | /** |
---|
185 | * At least as a workaround we store the irredundant facets of a matrix here. |
---|
186 | * This is needed to compute an interior points of a cone. Note that there |
---|
187 | * will be non-flippable facets in it! |
---|
188 | */ |
---|
189 | dd_MatrixPtr ddFacets; //Matrix to store irredundant facets of the cone |
---|
190 | |
---|
191 | /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/ |
---|
192 | ideal gcBasis; //GB of the cone, set by gcone::getGB(); |
---|
193 | gcone *next; //Pointer to next cone |
---|
194 | gcone *prev; |
---|
195 | |
---|
196 | gcone(); |
---|
197 | gcone(ring r, ideal I); |
---|
198 | gcone(const gcone& gc, const facet &f); |
---|
199 | ~gcone(); |
---|
200 | inline int getCounter(); |
---|
201 | inline ring getBaseRing(); |
---|
202 | inline void setIntPoint(intvec *iv); |
---|
203 | inline intvec *getIntPoint(bool shallow=FALSE); |
---|
204 | inline void showIntPoint(); |
---|
205 | inline void setNumFacets(); |
---|
206 | inline int getNumFacets(); |
---|
207 | inline int getUCN(); |
---|
208 | inline int getPredUCN(); |
---|
209 | volatile void showFacets(short codim=1); |
---|
210 | inline volatile void showSLA(facet &f); |
---|
211 | inline void idDebugPrint(const ideal &I); |
---|
212 | inline void invPrint(const ideal &I); |
---|
213 | inline bool isMonomial(const ideal &I); |
---|
214 | inline intvec *ivNeg(const intvec *iv); |
---|
215 | // inline int dotProduct(intvec &iva, intvec &ivb); |
---|
216 | inline int dotProduct(const intvec &iva, const intvec &ivb); |
---|
217 | inline bool isParallel(const intvec &a, const intvec &b); |
---|
218 | // inline int dotProduct(const intvec* a, const intvec *b); |
---|
219 | // inline bool isParallel(const intvec* a, const intvec* b); |
---|
220 | inline bool ivAreEqual(intvec &a, intvec &b); |
---|
221 | inline bool areEqual( facet *f, facet *g); |
---|
222 | inline bool areEqual2(facet* f, facet *g); |
---|
223 | inline int intgcd(const int &a, const int &b); |
---|
224 | inline void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE); |
---|
225 | inline void readConeFromFile(int gcNum, gcone *gc); |
---|
226 | inline intvec f2M(gcone *gc, facet *f, int n=1); |
---|
227 | inline void sortRays(gcone *gc); |
---|
228 | //The real stuff |
---|
229 | inline void getConeNormals(const ideal &I, bool compIntPoint=FALSE); |
---|
230 | inline void getCodim2Normals(const gcone &gc); |
---|
231 | inline void getExtremalRays(const gcone &gc); |
---|
232 | inline void flip(ideal gb, facet *f); |
---|
233 | inline void flip2(const ideal gb, facet *f); |
---|
234 | inline void computeInv(const ideal &gb, ideal &inv, const intvec &f);// poly restOfDiv(poly const &f, ideal const &I); removed with r12286 |
---|
235 | inline ideal ffG(const ideal &H, const ideal &G); |
---|
236 | inline void getGB(ideal const &inputIdeal); |
---|
237 | inline void interiorPoint( dd_MatrixPtr &M, intvec &iv); |
---|
238 | inline void interiorPoint2(); //removed Feb 8th, 2010, new method Feb 19th, 2010 |
---|
239 | inline void preprocessInequalities(dd_MatrixPtr &M); |
---|
240 | ring rCopyAndAddWeight(const ring &r, intvec *ivw); |
---|
241 | ring rCopyAndAddWeight2(const ring &, const intvec *, const intvec *); |
---|
242 | ring rCopyAndChangeWeight(const ring &r, intvec *ivw); |
---|
243 | // void reverseSearch(gcone *gcAct); //NOTE both removed from r12286 |
---|
244 | // bool isSearchFacet(gcone &gcTmp, facet *testfacet); |
---|
245 | void noRevS(gcone &gcRoot, bool usingIntPoint=FALSE); |
---|
246 | inline void makeInt(const dd_MatrixPtr &M, const int line, intvec &n); |
---|
247 | inline void normalize(); |
---|
248 | facet * enqueueNewFacets(facet *f); |
---|
249 | facet * enqueue2(facet *f); |
---|
250 | dd_MatrixPtr facets2Matrix(const gcone &gc); |
---|
251 | /** Compute the lineality space Ax=0 and return it as dd_MatrixPtr dd_LinealitySpace*/ |
---|
252 | inline dd_MatrixPtr computeLinealitySpace(); |
---|
253 | inline bool iv64isStrictlyPositive(const intvec *); |
---|
254 | /** Exchange 2 ordertype_a by just 1 */ |
---|
255 | inline void replaceDouble_ringorder_a_ByASingleOne(); |
---|
256 | // static void gcone::idPrint(ideal &I); |
---|
257 | friend class facet; |
---|
258 | }; |
---|
259 | lists lprepareResult(gcone *gc, const int n); |
---|
260 | // bool iv64isStrictlyPositive(intvec *); |
---|
261 | #endif |
---|
262 | #endif |
---|