source: git/kernel/gfan.h @ be7ab3f

spielwiese
Last change on this file since be7ab3f was be7ab3f, checked in by Martin Monerjan, 14 years ago
Some asserts Facet copy constructor: no ivCopy of f.interiorPoint, just referencing. Bad style, should go into some kind of shallow copy Skeleton for shallow copy getExtremalRays: For hom input, compute rays, add them up and add a suitable mulitple of (1,_1); for non hom we add the constraints for the pos orthant. Both ways we get strictly pos int points. flip2: Different approach to flipping using an ordering (omega,-N,dp); requires no marking and is appreciably faster than the old way interiorPoint2: Legacy code, no longer used since int points are now computed in getExtremalRays rCopyAndAddWeight2 (a(),a(),dp) Earlier computation of linspace and ext rays for 1st cone interiorPoint goes into the searchlist replaceDouble_ringorder_a_ByASingleOne() changes in readConeFromFile - to be undone git-svn-id: file:///usr/local/Singular/svn/trunk@12589 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 9.0 KB
Line 
1/*
2gfan.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
23extern int gfanHeuristic;
24// extern dd_MatrixPtr ddLinealitySpace;
25#define gfanp
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);
35lists gfan(ideal inputIdeal, int heuristic);
36//int dotProduct(intvec a, intvec b);
37//bool isParallel(intvec a, intvec b);
38
39class 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                ring flipRing;          //the ring on the other side of the facet
78                unsigned numRays;       //Number of spanning rays of the cone                   
79                /** The default constructor. */
80                facet();
81                /** Constructor for lower dimensional faces*/
82                facet(int const &n);
83                /**  The copy constructor */
84                facet(const facet& f);
85                facet(const facet& f, bool shallow);
86                /** The default destructor */
87                ~facet();
88                /** Comparison operator*/
89//              inline bool operator==(const facet *f,const facet *g);                 
90                /** \brief Comparison of facets*/
91                inline bool areEqual(facet *f, facet *g);
92                /** Stores the facet normal \param intvec*/
93                inline void setFacetNormal(intvec *iv);
94                /** Hopefully returns the facet normal */
95                inline intvec *getFacetNormal();
96                /** Return a reference to the facet normal*/
97                inline const intvec *getRef2FacetNormal();
98                /** Method to print the facet normal*/
99                inline void printNormal();
100                /** Store the flipped GB*/
101                inline void setFlipGB(ideal I);
102                /** Return the flipped GB*/
103                inline ideal getFlipGB();
104                /** Print the flipped GB*/
105                inline void printFlipGB();
106                /** Set the UCN */
107                inline void setUCN(int n);
108                /** \brief Get the UCN
109                 * Returns the UCN iff this != NULL, else -1
110                 */
111                inline int getUCN();
112                /** Store an interior point of the facet */
113                inline void setInteriorPoint(intvec *iv);
114                inline intvec *getInteriorPoint();
115                /** \brief Debugging function
116                 * prints the facet normal an all (codim-2)-facets that belong to it
117                 */
118                volatile void fDebugPrint();
119                friend class gcone;             
120};
121
122/**
123 *\brief Implements the cone structure
124 *
125 * A cone is represented by a linked list of facet normals
126 * @see facet
127 */
128
129class gcone
130{
131        private:               
132//              ring rootRing;          //good to know this -> generic walk
133                ideal inputIdeal;       //the original
134                ring baseRing;          //the basering of the cone                             
135                intvec *ivIntPt;        //an interior point of the cone
136                int UCN;                //unique number of the cone
137                int pred;               //UCN of the cone this one is derived from
138                static int counter;
139               
140        public: 
141                /** \brief Pointer to the first facet */
142                facet *facetPtr;        //Will hold the adress of the first facet; set by gcone::getConeNormals
143#ifdef gfanp
144                static float time_getConeNormals;
145                static float time_getCodim2Normals;
146                static float t_getExtremalRays;
147                static float time_flip;
148                static float time_flip2;
149                static float t_areEqual;
150                static float t_ffG;
151                static float t_markings;
152                static float t_dd;
153                static float t_kStd;
154                static float time_enqueue;             
155                static float time_computeInv;
156                static float t_ddMC;
157                static float t_mI;
158                static float t_iP;
159                static unsigned parallelButNotEqual;
160                static unsigned numberOfFacetChecks;
161#endif
162                /** Matrix to contain the homogeneity/lineality space */
163                static dd_MatrixPtr dd_LinealitySpace;
164                static int lengthOfSearchList;
165                /** Maximum size of the searchlist*/
166                static int maxSize;
167                /** is the ideal homogeneous? */
168                static bool hasHomInput;
169                /** # of variables in the ring */
170                int numVars;            //#of variables in the ring
171                /** The hilbert function - for the homogeneous case*/
172                static intvec *hilbertFunction;
173               
174                /** # of facets of the cone
175                 * This value is set by gcone::getConeNormals
176                 */
177                int numFacets;          //#of facets of the cone
178               
179                /**
180                 * At least as a workaround we store the irredundant facets of a matrix here.
181                 * This is needed to compute an interior points of a cone. Note that there
182                 * will be non-flippable facets in it!           
183                 */
184                dd_MatrixPtr ddFacets;  //Matrix to store irredundant facets of the cone
185               
186                /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
187                ideal gcBasis;          //GB of the cone, set by gcone::getGB();
188                gcone *next;            //Pointer to next cone
189                gcone *prev;
190               
191                gcone();
192                gcone(ring r, ideal I);
193                gcone(const gcone& gc, const facet &f);
194                ~gcone();
195                inline int getCounter();
196                inline ring getBaseRing();
197                inline void setIntPoint(intvec *iv);
198                inline intvec *getIntPoint();
199                inline void showIntPoint();
200                inline void setNumFacets();
201                inline int getNumFacets();
202                inline int getUCN();
203                inline int getPredUCN();               
204                volatile void showFacets(short codim=1);
205                inline volatile void showSLA(facet &f);
206                inline void idDebugPrint(const ideal &I);
207                inline void invPrint(const ideal &I);
208                inline bool isMonomial(const ideal &I);
209                inline intvec *ivNeg(const intvec *iv);
210//              inline int dotProduct(intvec &iva, intvec &ivb);
211                inline int dotProduct(const intvec &iva, const intvec &ivb);
212                inline bool isParallel(const intvec &a, const intvec &b);
213//              inline int dotProduct(const intvec* a, const intvec *b);
214//              inline bool isParallel(const intvec* a, const intvec* b);
215                inline bool areEqual(intvec &a, intvec &b);
216                inline bool areEqual( facet *f,  facet *g);
217                inline bool areEqual2(facet* f, facet *g);
218                inline int intgcd(int a, int b);
219                inline void writeConeToFile(const gcone &gc, bool usingIntPoints=FALSE);
220                inline void readConeFromFile(int gcNum, gcone *gc);
221                inline intvec f2M(gcone *gc, facet *f, int n=1);
222                inline void sortRays(gcone *gc);
223                //The real stuff
224                inline void getConeNormals(const ideal &I, bool compIntPoint=FALSE);
225                inline void getCodim2Normals(const gcone &gc);
226                inline void getExtremalRays(const gcone &gc);
227                inline void flip(ideal gb, facet *f);
228                inline void flip2(const ideal gb, facet *f);
229                inline void computeInv(const ideal &gb, ideal &inv, const intvec &f);//                 poly restOfDiv(poly const &f, ideal const &I); removed with r12286
230                inline ideal ffG(const ideal &H, const ideal &G);
231                inline void getGB(ideal const &inputIdeal);             
232                inline void interiorPoint( dd_MatrixPtr &M, intvec &iv);
233                inline void interiorPoint2(); //removed Feb 8th, 2010, new method Feb 19th, 2010
234                inline void preprocessInequalities(dd_MatrixPtr &M);
235                ring rCopyAndAddWeight(const ring &r, intvec *ivw);
236                ring rCopyAndAddWeight2(const ring &, const intvec *, const intvec *);
237                ring rCopyAndChangeWeight(const ring &r, intvec *ivw);         
238//              void reverseSearch(gcone *gcAct); //NOTE both removed from r12286
239//              bool isSearchFacet(gcone &gcTmp, facet *testfacet);
240                void noRevS(gcone &gcRoot, bool usingIntPoint=FALSE);
241                inline void makeInt(const dd_MatrixPtr &M, const int line, intvec &n);
242                inline void normalize();
243                facet * enqueueNewFacets(facet *f);
244                facet * enqueue2(facet *f);
245                dd_MatrixPtr facets2Matrix(const gcone &gc);
246                /** Compute the lineality space Ax=0 and return it as dd_MatrixPtr dd_LinealitySpace*/
247                inline dd_MatrixPtr computeLinealitySpace();
248                inline bool iv64isStrictlyPositive(const intvec *);
249                /** Exchange 2 ordertype_a by just 1 */
250                inline void replaceDouble_ringorder_a_ByASingleOne();
251//              static void gcone::idPrint(ideal &I);           
252                friend class facet;     
253};
254lists lprepareResult(gcone *gc, int n);
255// bool iv64isStrictlyPositive(intvec *);
256#endif
257#endif
Note: See TracBrowser for help on using the repository browser.