source: git/kernel/gfan.cc @ 5435c3

spielwiese
Last change on this file since 5435c3 was 5435c3, checked in by Martin Monerjan, 15 years ago
UCN Setters and getters method gcone:noRevS git-svn-id: file:///usr/local/Singular/svn/trunk@11792 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 43.5 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
4$Date: 2009-05-12 15:30:29 $
5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.50 2009-05-12 15:30:29 monerjan Exp $
6$Id: gfan.cc,v 1.50 2009-05-12 15:30:29 monerjan Exp $
7*/
8
9#include "mod2.h"
10
11#ifdef HAVE_GFAN
12
13#include "kstd1.h"
14#include "kutil.h"
15#include "intvec.h"
16#include "polys.h"
17#include "ideals.h"
18#include "kmatrix.h"
19#include "fast_maps.h"  //Mapping of ideals
20#include "maps.h"
21#include "ring.h"
22#include "prCopy.h"
23#include <iostream>
24#include <bitset>
25#include <fstream>      //read-write cones to files
26
27/*DO NOT REMOVE THIS*/
28#ifndef GMPRATIONAL
29#define GMPRATIONAL
30#endif
31
32//Hacks for different working places
33#define ITWM
34
35#ifdef UNI
36#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
37#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
38#endif
39
40#ifdef HOME
41#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
42#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
43#endif
44
45#ifdef ITWM
46#include "/u/slg/monerjan/cddlib/include/setoper.h"
47#include "/u/slg/monerjan/cddlib/include/cdd.h"
48#include "/u/slg/monerjan/cddlib/include/cddmp.h"
49#endif
50
51#ifndef gfan_DEBUG
52#define gfan_DEBUG
53#endif
54
55//#include gcone.h
56using namespace std;
57
58/**
59*\brief Class facet
60*       Implements the facet structure as a linked list
61*
62*/
63class facet
64{
65        private:
66                /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
67                intvec *fNormal;
68               
69                /** \brief An interior point of the facet*/
70                intvec *interiorPoint;
71               
72                /** \brief Universal Cone Number
73                * The number of the cone the facet belongs to
74                */
75                int UCN;
76               
77                /** \brief The Groebner basis on the other side of a shared facet
78                 *
79                 * In order not to have to compute the flipped GB twice we store the basis we already get
80                 * when identifying search facets. Thus in the next step of the reverse search we can
81                 * just copy the old cone and update the facet and the gcBasis.
82                 * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
83                 */
84                ideal flipGB;           //The Groebner Basis on the other side, computed via gcone::flip
85
86                       
87        public:         
88                //bool isFlippable;     //flippable facet? Want to have cone->isflippable.facet[i]
89                bool isIncoming;        //Is the facet incoming or outgoing?
90                facet *next;            //Pointer to next facet
91               
92                /** The default constructor. Do I need a constructor of type facet(intvec)? */
93                facet()                 
94                {
95                        // Pointer to next facet.  */
96                        /* Defaults to NULL. This way there is no need to check explicitly */
97                        this->next=NULL; 
98                }
99               
100                /** The default destructor */
101                ~facet(){;}
102               
103                /** Stores the facet normal \param intvec*/
104                void setFacetNormal(intvec *iv)
105                {
106                        this->fNormal = ivCopy(iv);
107                        //return;
108                }
109               
110                /** Hopefully returns the facet normal */
111                intvec *getFacetNormal()
112                {       
113                        //this->fNormal->show();        cout << endl;   
114                        return this->fNormal;
115                }
116               
117                /** Method to print the facet normal*/
118                void printNormal()
119                {
120                        fNormal->show();
121                }
122               
123                /** Store the flipped GB*/
124                void setFlipGB(ideal I)
125                {
126                        this->flipGB=I;
127                }
128               
129                /** Return the flipped GB*/
130                ideal getFlipGB()
131                {
132                        return this->flipGB;
133                }
134               
135                /** Print the flipped GB*/
136                void printFlipGB()
137                {
138                        idShow(this->flipGB);
139                }
140               
141                void setUCN(int n)
142                {
143                        this->UCN=n;
144                }
145               
146                int getUCN()
147                {
148                        return this->UCN;
149                }
150               
151                void setInteriorPoint(intvec *iv)
152                {
153                        this->interiorPoint = ivCopy(iv);
154                }
155               
156                intvec *getInteriorPoint()
157                {
158                        return this->interiorPoint;
159                }
160               
161                /*bool isFlippable(intvec &load)
162                {
163                        bool res=TRUE;                 
164                        int jj;
165                        for (int jj = 0; jj<load.length(); jj++)
166                        {
167                                intvec *ivCanonical = new intvec(load.length());
168                                (*ivCanonical)[jj]=1;                           
169                                if (ivMult(&load,ivCanonical)<0)
170                                {
171                                        res=FALSE;
172                                        break;
173                                }
174                        }
175                        return res;
176                       
177                        /*while (dotProduct(load,ivCanonical)>=0)
178                        {
179                                if (jj!=this->numVars)
180                                {
181                                        intvec *ivCanonical = new intvec(this->numVars);
182                                        (*ivCanonical)[jj]=1;                   
183                                        res=TRUE;
184                                        jj += 1;
185                                }
186                        }
187                        if (jj==this->numVars)
188                        {                       
189                                delete ivCanonical;
190                                return FALSE;
191                        }
192                        else
193                        {
194                                delete ivCanonical;
195                                return TRUE;
196                }*/                                             
197                //}//bool isFlippable(facet &f)
198               
199               
200                friend class gcone;     //Bad style
201};
202
203/**
204*\brief Implements the cone structure
205*
206* A cone is represented by a linked list of facet normals
207* @see facet
208*/
209/*class gcone
210finally this should become s.th. like gconelib.{h,cc} to provide an API
211*/
212class gcone
213{
214        private:
215                int numFacets;          //#of facets of the cone
216                ring rootRing;          //good to know this -> generic walk
217                ideal inputIdeal;       //the original
218                ring baseRing;          //the basering of the cone             
219                /* TODO in order to save memory use pointers to rootRing and inputIdeal instead */
220                intvec *ivIntPt;        //an interior point of the cone
221                int UCN;                //unique number of the cone
222               
223        public: 
224                /** Pointer to the first facet */
225                facet *facetPtr;        //Will hold the adress of the first facet; set by gcone::getConeNormals
226               
227                /** # of variables in the ring */
228                int numVars;            //#of variables in the ring
229               
230                /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
231                ideal gcBasis;          //GB of the cone, set by gcone::getGB();
232                gcone *next;            //Pointer to *previous* cone in search tree     
233               
234                /** \brief Default constructor.
235                *
236                * Initialises this->next=NULL and this->facetPtr=NULL
237                */
238                gcone()
239                {
240                        this->next=NULL;
241                        this->facetPtr=NULL;
242                        this->baseRing=currRing;
243                        this->UCN=1;
244                }
245               
246                /** \brief Constructor with ring and ideal
247                *
248                * This constructor takes the root ring and the root ideal as parameters and stores
249                * them in the private members gcone::rootRing and gcone::inputIdeal
250                * Since knowledge of the root ring is only needed when using reverse search,
251                * this constructor is not needed when using the "second" method
252                */
253                gcone(ring r, ideal I)
254                {
255                        this->next=NULL;
256                        this->facetPtr=NULL;
257                        this->rootRing=r;
258                        this->inputIdeal=I;
259                        this->baseRing=currRing;
260                        this->UCN=1;
261                }
262               
263                /** \brief Copy constructor
264                *
265                * Copies one cone, sets this->gcBasis to the flipped GB and reverses the
266                * direction of the according facet normal
267                */                     
268                gcone(const gcone& gc)         
269                {
270                        this->next=NULL;
271                        this->numVars=gc.numVars;
272                        this->UCN=(gc.UCN)+1;   //add 1 to the UCN of previous cone
273                        facet *fAct= new facet();                       
274                        this->facetPtr=fAct;
275                       
276                        intvec *ivtmp = new intvec(this->numVars);                                             
277                        ivtmp = gc.facetPtr->getFacetNormal();
278                        //ivtmp->show();                       
279                       
280                        ideal gb;
281                        gb=gc.facetPtr->getFlipGB();                   
282                        this->gcBasis=gb;//gc.facetPtr->getFlipGB();    //this cone's GB is the flipped GB
283                        //idShow(this->gcBasis);
284                       
285                        /*Reverse direction of the facet normal to make it an inner normal*/                   
286                        for (int ii=0; ii<this->numVars;ii++)
287                        {
288                                (*ivtmp)[ii]=-(*ivtmp)[ii];                             
289                        }
290                       
291                        fAct->setFacetNormal(ivtmp);
292                        delete ivtmp;                                           
293                }
294               
295                /** \brief Default destructor */
296                ~gcone(){;}             //destructor
297                       
298               
299                /** \brief Set the interior point of a cone */
300                void setIntPoint(intvec *iv)
301                {
302                        this->ivIntPt=ivCopy(iv);
303                }
304               
305                /** \brief Return the interior point */
306                intvec *getIntPoint()
307                {
308                        return this->ivIntPt;
309                }
310               
311                void showIntPoint()
312                {
313                        ivIntPt->show();
314                }
315               
316                /** \brief Set gcone::numFacets */
317                void setNumFacets()
318                {
319                }
320               
321                /** \brief Get gcone::numFacets */
322                int getNumFacets()
323                {
324                        return this->numFacets;
325                }
326               
327                int getUCN()
328                {
329                        return this->UCN;
330                }
331               
332                /** \brief Compute the normals of the cone
333                *
334                * This method computes a representation of the cone in terms of facet normals. It takes an ideal
335                * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
336                * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
337                * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
338                * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
339                * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
340                * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
341                *
342                * Optionally, if the parameter bool compIntPoint is set to TRUE the method will also compute
343                * an interior point of the cone.
344                */
345                void getConeNormals(ideal const &I, bool compIntPoint=FALSE)
346                {
347#ifdef gfan_DEBUG
348                        std::cout << "*** Computing Inequalities... ***" << std::endl;
349#endif         
350                        //All variables go here - except ineq matrix and *v, see below
351                        int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
352                        int pCompCount;                 // # of terms in a poly
353                        poly aktpoly;
354                        int numvar = pVariables;        // # of variables in a polynomial (or ring?)
355                        int leadexp[numvar];            // dirty hack of exp.vects
356                        int aktexp[numvar];
357                        int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
358                        dd_rowrange ddrows;
359                        dd_colrange ddcols;
360                        dd_rowset ddredrows;            // # of redundant rows in ddineq
361                        dd_rowset ddlinset;             // the opposite
362                        dd_rowindex ddnewpos;             // all to make dd_Canonicalize happy
363                        dd_NumberType ddnumb=dd_Integer; //Number type
364                        dd_ErrorType dderr=dd_NoError;  //
365                        // End of var declaration
366#ifdef gfan_DEBUG
367                        cout << "The Groebner basis has " << lengthGB << " elements" << endl;
368                        cout << "The current ring has " << numvar << " variables" << endl;
369#endif
370                        cols = numvar;
371               
372                        //Compute the # inequalities i.e. rows of the matrix
373                        rows=0; //Initialization
374                        for (int ii=0;ii<IDELEMS(I);ii++)
375                        {
376                                aktpoly=(poly)I->m[ii];
377                                rows=rows+pLength(aktpoly)-1;
378                        }
379#ifdef gfan_DEBUG
380                        cout << "rows=" << rows << endl;
381                        cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
382#endif
383                        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
384                        dd_set_global_constants();
385                        ddrows=rows;
386                        ddcols=cols;
387                        dd_MatrixPtr ddineq;            //Matrix to store the inequalities                     
388                        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
389               
390                        // We loop through each g\in GB and compute the resulting inequalities
391                        for (int i=0; i<IDELEMS(I); i++)
392                        {
393                                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
394                                pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
395                                cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
396               
397                                int *v=(int *)omAlloc((numvar+1)*sizeof(int));
398                                pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
399                               
400                                //Store leadexp for aktpoly
401                                for (int kk=0;kk<numvar;kk++)
402                                {
403                                        leadexp[kk]=v[kk+1];
404                                        //Since we need to know the difference of leadexp with the other expvects we do nothing here
405                                        //but compute the diff below
406                                }
407               
408                               
409                                while (pNext(aktpoly)!=NULL) //move to next term until NULL
410                                {
411                                        aktpoly=pNext(aktpoly);
412                                        pSetm(aktpoly);         //doesn't seem to help anything
413                                        pGetExpV(aktpoly,v);
414                                        for (int kk=0;kk<numvar;kk++)
415                                        {
416                                                aktexp[kk]=v[kk+1];
417                                                //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
418                                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
419                                        }
420                                        aktmatrixrow=aktmatrixrow+1;
421                                }//while
422               
423                        } //for
424               
425                        //Maybe add another row to contain the constraints of the standard simplex?
426
427#ifdef gfan_DEBUG
428                        cout << "The inequality matrix is" << endl;
429                        dd_WriteMatrix(stdout, ddineq);
430#endif
431
432                        // The inequalities are now stored in ddineq
433                        // Next we check for superflous rows
434                        ddredrows = dd_RedundantRows(ddineq, &dderr);
435                        if (dderr!=dd_NoError)                  // did an error occur?
436                        {
437                                dd_WriteErrorMessages(stderr,dderr);    //if so tell us
438                        } else
439                        {
440                                cout << "Redundant rows: ";
441                                set_fwrite(stdout, ddredrows);          //otherwise print the redundant rows
442                        }//if dd_Error
443               
444                        //Remove reduntant rows here!
445                        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
446                        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
447                        ddcols = ddineq->colsize;
448#ifdef gfan_DEBUG
449                        cout << "Having removed redundancies, the normals now read:" << endl;
450                        dd_WriteMatrix(stdout,ddineq);
451                        cout << "Rows = " << ddrows << endl;
452                        cout << "Cols = " << ddcols << endl;
453#endif
454                       
455                        /*Write the normals into class facet*/
456#ifdef gfan_DEBUG
457                        cout << "Creating list of normals" << endl;
458#endif
459                        /*The pointer *fRoot should be the return value of this function*/
460                        facet *fRoot = new facet();     //instantiate new facet
461                        this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
462                        facet *fAct;                    //instantiate pointer to active facet
463                        fAct = fRoot;                   //Seems to do the trick. fRoot and fAct have to point to the same adress!
464                        //std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
465                        for (int kk = 0; kk<ddrows; kk++)
466                        {
467                                intvec *load = new intvec(this->numVars);       //intvec to store a single facet normal that will then be stored via setFacetNormal
468                                for (int jj = 1; jj <ddcols; jj++)
469                                {
470#ifdef GMPRATIONAL
471                                        double foo;
472                                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
473#ifdef gfan_DEBUG
474                                        std::cout << "fAct is " << foo << " at " << fAct << std::endl;
475#endif
476                                        (*load)[jj-1] = (int)foo;               //store typecasted entry at pos jj-1 of load
477#endif                 
478
479                                        //(*load)[jj-1] = (int)foo;             //store typecasted entry at pos jj-1 of load
480                                        //check for flipability here
481//                                      if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
482//                                      {
483                                                //fAct->next = new facet();     //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                               
484//                                      }
485                                }//for (int jj = 1; jj <ddcols; jj++)
486                                /*Quick'n'dirty hack for flippability*/ 
487                                bool isFlippable=FALSE;                         
488                                for (int jj = 0; jj<load->length(); jj++)
489                                {
490                                        intvec *ivCanonical = new intvec(load->length());
491                                        (*ivCanonical)[jj]=1;
492                                        cout << "dotProd=" << dotProduct(*load,*ivCanonical) << endl;
493                                        if (dotProduct(*load,*ivCanonical)<0)                                   
494                                        {
495                                                isFlippable=TRUE;
496                                                break;  //URGHS
497                                        }
498                                }       
499                                if (isFlippable==FALSE)
500                                {
501                                        cout << "Ignoring facet";
502                                        load->show();
503                                        //fAct->next=NULL;
504                                }
505                                else
506                                {       /*Now load should be full and we can call setFacetNormal*/
507                                        fAct->setFacetNormal(load);
508                                        fAct->next = new facet();                                       
509                                        fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
510                                }//if (isFlippable==FALSE)
511                                delete load;
512                        }//for (int kk = 0; kk<ddrows; kk++)
513                        /*
514                        Now we should have a linked list containing the facet normals of those facets that are
515                        -irredundant
516                        -flipable
517                        Adressing is done via *facetPtr
518                        */
519                       
520                        if (compIntPoint==TRUE)
521                        {
522                                intvec *iv = new intvec(this->numVars);
523                                interiorPoint(ddineq, *iv);     //NOTE ddineq contains non-flippable facets
524                                this->setIntPoint(iv);  //stores the interior point in gcone::ivIntPt
525                                delete iv;
526                        }
527                       
528                        //Clean up but don't delete the return value! (Whatever it will turn out to be)                 
529                        dd_FreeMatrix(ddineq);
530                        set_free(ddredrows);
531                        free(ddnewpos);
532                        set_free(ddlinset);
533                        dd_free_global_constants();
534
535                }//method getConeNormals(ideal I)       
536               
537               
538                /** \brief Compute the Groebner Basis on the other side of a shared facet
539                *
540                * Implements algorithm 4.3.2 from Anders' thesis.
541                * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
542                * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
543                * is parallel to \f$ leadexp(g) \f$
544                * Parallelity is checked using basic linear algebra. See gcone::isParallel.
545                * Other possibilities include computing the rank of the matrix consisting of the vectors in question and
546                * computing an interior point of the facet and taking all terms having the same weight with respect
547                * to this interior point.
548                *\param ideal, facet
549                * Input: a marked,reduced Groebner basis and a facet
550                */
551                void flip(ideal gb, facet *f)           //Compute "the other side"
552                {                       
553                        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
554                        fNormal = f->getFacetNormal();  //read this->fNormal;
555#ifdef gfan_DEBUG
556                        std::cout << "===" << std::endl;
557                        std::cout << "running gcone::flip" << std::endl;
558//                      std::cout << "fNormal=";
559//                      fNormal->show();
560//                      std::cout << std::endl;
561#endif                         
562                        /*1st step: Compute the initial ideal*/
563                        poly initialFormElement[IDELEMS(gb)];   //array of #polys in GB to store initial form
564                        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
565                        poly aktpoly;
566                        intvec *check = new intvec(this->numVars);      //array to store the difference of LE and v
567                       
568                        for (int ii=0;ii<IDELEMS(gb);ii++)
569                        {
570                                aktpoly = (poly)gb->m[ii];                                                             
571                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
572                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
573                                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
574                                initialFormElement[ii]=pHead(aktpoly);
575                               
576                                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
577                                {
578                                        aktpoly=pNext(aktpoly); //next term
579                                        pSetm(aktpoly);
580                                        pGetExpV(aktpoly,v);           
581                                        /* Convert (int)v into (intvec)check */                 
582                                        for (int jj=0;jj<this->numVars;jj++)
583                                        {
584                                                //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
585                                                //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
586                                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
587                                        }
588#ifdef gfan_DEBUG
589//                                      cout << "check=";                       
590//                                      check->show();
591//                                      cout << endl;
592#endif
593                                        //TODO why not *check, *fNormal????
594                                        if (isParallel(*check,*fNormal)) //pass *check when
595                                        {
596//                                              cout << "Parallel vector found, adding to initialFormElement" << endl;                 
597                                                initialFormElement[ii] = pAdd(pCopy(initialFormElement[ii]),(poly)pHead(aktpoly));
598                                        }                                               
599                                }//while
600#ifdef gfan_DEBUG
601                                cout << "Initial Form=";                               
602                                pWrite(initialFormElement[ii]);
603                                cout << "---" << endl;
604#endif
605                                /*Now initialFormElement must be added to (ideal)initialForm */
606                                initialForm->m[ii]=initialFormElement[ii];
607                        }//for                 
608#ifdef gfan_DEBUG
609                        cout << "Initial ideal is: " << endl;
610                        idShow(initialForm);
611                        //f->printFlipGB();
612                        cout << "===" << endl;
613#endif
614                        //delete check;
615                       
616                        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
617                        /*Substep 2.1
618                        compute $G_{-\alpha}(in_v(I))
619                        see journal p. 66
620                        */
621                        ring srcRing=currRing;
622
623                        //intvec *negfNormal = new intvec(this->numVars);
624                        //negfNormal=ivNeg(fNormal);
625                        ring tmpRing=rCopyAndAddWeight(srcRing,ivNeg(fNormal));
626                        rChangeCurrRing(tmpRing);
627                       
628                        //rWrite(currRing); cout << endl;
629                       
630                        ideal ina;                     
631                        ina=idrCopyR(initialForm,srcRing);                     
632#ifdef gfan_DEBUG
633                        cout << "ina=";
634                        idShow(ina); cout << endl;
635#endif
636                        ideal H;
637                        H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
638                        idSkipZeroes(H);
639#ifdef gfan_DEBUG
640//                      cout << "H="; idShow(H); cout << endl;
641#endif
642                        /*Substep 2.2
643                        do the lifting and mark according to H
644                        */
645                        rChangeCurrRing(srcRing);
646                        ideal srcRing_H;
647                        ideal srcRing_HH;                       
648                        srcRing_H=idrCopyR(H,tmpRing);
649#ifdef gfan_DEBUG
650//                      cout << "srcRing_H = ";
651//                      idShow(srcRing_H); cout << endl;
652#endif
653                        srcRing_HH=ffG(srcRing_H,this->gcBasis);               
654#ifdef gfan_DEBUG
655//                      cout << "srcRing_HH = ";
656//                      idShow(srcRing_HH); cout << endl;
657#endif
658                        /*Substep 2.2.1
659                        Mark according to G_-\alpha
660                        Here we have a minimal basis srcRing_HH. In order to turn this basis into a reduced basis
661                        we have to compute an interior point of C(srcRing_HH). For this we need to know the cone
662                        represented by srcRing_HH MARKED ACCORDING TO G_{-\alpha}
663                        Thus we check whether the leading monomials of srcRing_HH and srcRing_H coincide. If not we
664                        compute the difference accordingly
665                        */
666                        dd_set_global_constants();
667                        bool markingsAreCorrect=FALSE;
668                        dd_MatrixPtr intPointMatrix;
669                        int iPMatrixRows=0;
670                        dd_rowrange aktrow=0;                   
671                        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
672                        {
673                                poly aktpoly=(poly)srcRing_HH->m[ii];
674                                iPMatrixRows = iPMatrixRows+pLength(aktpoly)-1;
675                        }
676                        /* additionally one row for the standard-simplex and another for a row that becomes 0 during
677                        construction of the differences
678                        */
679                        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
680                        intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
681                       
682                        for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
683                        {
684                                markingsAreCorrect=FALSE;       //crucial to initialise here
685                                poly aktpoly=srcRing_HH->m[ii];
686                                /*Comparison of leading monomials is done via exponent vectors*/
687                                for (int jj=0;jj<IDELEMS(H);jj++)
688                                {
689                                        int *src_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
690                                        int *dst_ExpV = (int *)omAlloc((this->numVars+1)*sizeof(int));
691                                        pGetExpV(aktpoly,src_ExpV);
692                                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
693                                        pGetExpV(pCopy(H->m[ii]),dst_ExpV);
694                                        rChangeCurrRing(srcRing);
695                                        bool expVAreEqual=TRUE;
696                                        for (int kk=1;kk<=this->numVars;kk++)
697                                        {
698#ifdef gfan_DEBUG
699                                                //cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
700#endif
701                                                if (src_ExpV[kk]!=dst_ExpV[kk])
702                                                {
703                                                        expVAreEqual=FALSE;
704                                                }
705                                        }                                       
706                                        //if (*src_ExpV == *dst_ExpV)
707                                        if (expVAreEqual==TRUE)
708                                        {
709                                                markingsAreCorrect=TRUE; //everything is fine
710#ifdef gfan_DEBUG
711//                                              cout << "correct markings" << endl;
712#endif
713                                        }//if (pHead(aktpoly)==pHead(H->m[jj])
714                                        delete src_ExpV;
715                                        delete dst_ExpV;
716                                }//for (int jj=0;jj<IDELEMS(H);jj++)
717                               
718                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
719                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
720                                if (markingsAreCorrect==TRUE)
721                                {
722                                        pGetExpV(aktpoly,leadExpV);
723                                }
724                                else
725                                {
726                                        rChangeCurrRing(tmpRing);
727                                        pGetExpV(pHead(H->m[ii]),leadExpV); //We use H->m[ii] as leading monomial
728                                        rChangeCurrRing(srcRing);
729                                }
730                                /*compute differences of the expvects*/                         
731                                while (pNext(aktpoly)!=NULL)
732                                {
733                                        /*The following if-else-block makes sure the first term (i.e. the wrongly marked term)
734                                        is not omitted when computing the differences*/
735                                        if(markingsAreCorrect==TRUE)
736                                        {
737                                                aktpoly=pNext(aktpoly);
738                                                pGetExpV(aktpoly,v);
739                                        }
740                                        else
741                                        {
742                                                pGetExpV(pHead(aktpoly),v);
743                                                markingsAreCorrect=TRUE;
744                                        }
745
746                                        for (int jj=0;jj<this->numVars;jj++)
747                                        {                               
748                                                /*Store into ddMatrix*/                                                         
749                                                dd_set_si(intPointMatrix->matrix[aktrow][jj+1],leadExpV[jj+1]-v[jj+1]);
750                                        }
751                                        aktrow +=1;
752                                }
753                                delete v;
754                                delete leadExpV;
755                        }//for (int ii=0;ii<IDELEMS(srcRing_HH);ii++)
756                        /*Now we add the constraint for the standard simplex*/
757                        /*NOTE:Might actually work without the standard simplex*/
758                        dd_set_si(intPointMatrix->matrix[aktrow][0],-1);
759                        for (int jj=1;jj<=this->numVars;jj++)
760                        {
761                                dd_set_si(intPointMatrix->matrix[aktrow][jj],1);
762                        }
763                        dd_WriteMatrix(stdout,intPointMatrix);
764                        intvec *iv_weight = new intvec(this->numVars);
765                        interiorPoint(intPointMatrix, *iv_weight);      //iv_weight now contains the interior point
766                        dd_FreeMatrix(intPointMatrix);
767                        dd_free_global_constants();
768                       
769                        /*Step 3
770                        turn the minimal basis into a reduced one
771                        */
772                        int i,j;
773                        ring dstRing=rCopy0(srcRing);
774                        i=rBlocks(srcRing);
775                       
776                        dstRing->order=(int *)omAlloc((i+1)*sizeof(int));
777                        for(j=i;j>0;j--)
778                        {
779                                dstRing->order[j]=srcRing->order[j-1];
780                                dstRing->block0[j]=srcRing->block0[j-1];
781                                dstRing->block1[j]=srcRing->block1[j-1];
782                                if (srcRing->wvhdl[j-1] != NULL)
783                                {
784                                        dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
785                                }
786                        }
787                        dstRing->order[0]=ringorder_a;
788                        dstRing->order[1]=ringorder_dp;
789                        dstRing->order[2]=ringorder_C;                 
790                        dstRing->wvhdl[0] =( int *)omAlloc((iv_weight->length())*sizeof(int));
791                       
792                        for (int ii=0;ii<this->numVars;ii++)
793                        {                               
794                                dstRing->wvhdl[0][ii]=(*iv_weight)[ii];                         
795                        }
796                        rComplete(dstRing);
797                       
798                        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
799                        // Thus:
800                        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
801                        //cout << "PLING" << endl;
802                        /*ring dstRing=rCopy0(srcRing);
803                        for (int ii=0;ii<this->numVars;ii++)
804                        {                               
805                                dstRing->wvhdl[0][ii]=(*iv_weight)[ii];                         
806                        }*/
807                        rChangeCurrRing(dstRing);
808#ifdef gfan_DEBUG
809                        rWrite(dstRing); cout << endl;
810#endif
811                        ideal dstRing_I;                       
812                        dstRing_I=idrCopyR(srcRing_HH,srcRing);                 
813                        //validOpts<1>=TRUE;
814#ifdef gfan_DEBUG
815                        //idShow(dstRing_I);
816#endif                 
817                        BITSET save=test;
818                        test|=Sy_bit(OPT_REDSB);
819                        test|=Sy_bit(6);        //OPT_DEBUG                                     
820                        dstRing_I=kStd(idrCopyR(this->inputIdeal,this->rootRing),NULL,testHomog,NULL);                                 
821                        kInterRed(dstRing_I);
822                        idSkipZeroes(dstRing_I);
823                        test=save;
824                        /*End of step 3 - reduction*/
825                       
826                        f->setFlipGB(dstRing_I);//store the flipped GB
827#ifdef gfan_DEBUG
828                        cout << "Flipped GB is: " << endl;
829                        f->printFlipGB();
830#endif                 
831                }//void flip(ideal gb, facet *f)
832                               
833                /** \brief Compute the remainder of a polynomial by a given ideal
834                *
835                * Compute \f$ f^{\mathcal{G}} \f$
836                * Algorithm is taken from Cox, Little, O'Shea, IVA 2nd Ed. p 62
837                * However, since we are only interested in the remainder, there is no need to
838                * compute the factors \f$ a_i \f$
839                */
840                //NOTE: Should be replaced by kNF or kNF2
841                poly restOfDiv(poly const &f, ideal const &I)
842                {
843                        cout << "Entering restOfDiv" << endl;
844                        poly p=f;
845                        //pWrite(p);
846                        //poly r=kCreateZeroPoly(,currRing,currRing);   //The 0-polynomial, hopefully
847                        poly r=NULL;    //The zero polynomial
848                        int ii;
849                        bool divOccured;
850                       
851                        while (p!=NULL)
852                        {
853                                ii=1;
854                                divOccured=FALSE;
855                               
856                                while( (ii<=IDELEMS(I) && (divOccured==FALSE) ))
857                                {                                       
858                                        if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
859                                        {                                               
860                                                poly step1,step2,step3;
861                                                //cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
862                                                step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
863                                                //cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;                               
864                                                step2 = ppMult_qq(step1, I->m[ii-1]);                                           
865                                                step3 = pSub(pCopy(p), step2);
866                                                //p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));                       
867                                                //pSetm(p);
868                                                pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
869                                                p=step3;
870                                                //pWrite(p);                                           
871                                                divOccured=TRUE;
872                                        }
873                                        else
874                                        {
875                                                ii += 1;
876                                        }//if (pLmDivisibleBy(I->m[ii],p,currRing))
877                                }//while( (ii<IDELEMS(I) && (divOccured==FALSE) ))
878                                if (divOccured==FALSE)
879                                {
880                                        //cout << "TICK 5" << endl;
881                                        r=pAdd(pCopy(r),pHead(p));
882                                        pSetm(r);
883                                        pSort(r);
884                                        //cout << "r="; pWrite(r); cout << endl;
885                                       
886                                        if (pLength(p)!=1)
887                                        {
888                                                p=pSub(pCopy(p),pHead(p));      //Here it may occur that p=0 instead of p=NULL
889                                        }
890                                        else
891                                        {
892                                                p=NULL; //Hack to correct this situation                                               
893                                        }                                       
894                                        //cout << "p="; pWrite(p);
895                                }//if (divOccured==FALSE)
896                        }//while (p!=0)
897                        return r;
898                }//poly restOfDiv(poly const &f, ideal const &I)
899               
900                /** \brief Compute \f$ f-f^{\mathcal{G}} \f$
901                */
902                //NOTE: use kNF or kNF2 instead of restOfDivision
903                ideal ffG(ideal const &H, ideal const &G)
904                {
905                        cout << "Entering ffG" << endl;
906                        int size=IDELEMS(H);
907                        ideal res=idInit(size,1);
908                        poly temp1, temp2, temp3;       //polys to temporarily store values for pSub
909                        for (int ii=0;ii<size;ii++)
910                        {
911                                res->m[ii]=restOfDiv(H->m[ii],G);
912                                //res->m[ii]=kNF(H->m[ii],G);
913                                temp1=H->m[ii];
914                                temp2=res->m[ii];                               
915                                temp3=pSub(temp1, temp2);
916                                res->m[ii]=temp3;
917                                //res->m[ii]=pSub(temp1,temp2); //buggy
918                                //pSort(res->m[ii]);
919                                //pSetm(res->m[ii]);
920                                //cout << "res->m["<<ii<<"]=";pWrite(res->m[ii]);                                               
921                        }                       
922                        return res;
923                }
924               
925                /** \brief Compute a Groebner Basis
926                *
927                * Computes the Groebner basis and stores the result in gcone::gcBasis
928                *\param ideal
929                *\return void
930                */
931                void getGB(ideal const &inputIdeal)
932                {
933                        ideal gb;
934                        gb=kStd(inputIdeal,NULL,testHomog,NULL);
935                        idSkipZeroes(gb);
936                        this->gcBasis=gb;       //write the GB into gcBasis
937                }//void getGB
938               
939                /** \brief The Generic Groebner Walk due to FJLT
940                * Needed for computing the search facet
941                */
942                ideal GenGrbWlk(ideal, ideal)
943                {
944                }//GGW
945               
946                /** \brief Compute the negative of a given intvec
947                */             
948                intvec *ivNeg(intvec const &iv)
949                {
950                        intvec *res = new intvec(this->numVars);
951                        for(int ii=0;ii<this->numVars;ii++)
952                        {
953                                res[ii] = -(int)iv[ii];
954                        }
955                        return res;
956                }
957
958
959                /** \brief Compute the dot product of two intvecs
960                *
961                */
962                int dotProduct(intvec const &iva, intvec const &ivb)                           
963                {
964                        //intvec iva=a;
965                        //intvec ivb=b;
966                        int res=0;
967                        for (int i=0;i<this->numVars;i++)
968                        {
969                                res = res+(iva[i]*ivb[i]);
970                        }
971                        return res;
972                }//int dotProduct
973
974                /** \brief Check whether two intvecs are parallel
975                *
976                * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
977                */
978                bool isParallel(intvec const &a, intvec const &b)
979                {                       
980                        int lhs,rhs;
981                        lhs=dotProduct(a,b)*dotProduct(a,b);
982                        rhs=dotProduct(a,a)*dotProduct(b,b);
983                        //cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
984                        if (lhs==rhs)
985                        {
986                                return TRUE;
987                        }
988                        else
989                        {
990                                return FALSE;
991                        }
992                }//bool isParallel
993               
994                /** \brief Compute an interior point of a given cone
995                * Result will be written into intvec iv
996                */
997                void interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
998                {
999                        dd_LPPtr lp,lpInt;
1000                        dd_ErrorType err=dd_NoError;
1001                        dd_LPSolverType solver=dd_DualSimplex;
1002                        dd_LPSolutionPtr lpSol=NULL;
1003                        dd_rowset ddlinset,ddredrows;   //needed for dd_FindRelativeInterior
1004                        dd_rowindex ddnewpos;
1005                        dd_NumberType numb;     
1006                        //M->representation=dd_Inequality;
1007                        //M->objective-dd_LPMin;  //Not sure whether this is needed
1008                       
1009                        //NOTE: Make this n-dimensional!
1010                        //dd_set_si(M->rowvec[0],1);dd_set_si(M->rowvec[1],1);dd_set_si(M->rowvec[2],1);
1011                                                                       
1012                        //dd_MatrixCanonicalize(&M, &ddlinset, &ddredrows, &ddnewpos, &err);
1013                        //if (err!=dd_NoError){cout << "Error during dd_MatrixCanonicalize" << endl;}
1014                        //cout << "Tick 2" << endl;
1015                        //dd_WriteMatrix(stdout,M);
1016                       
1017                        lp=dd_Matrix2LP(M, &err);
1018                        if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}                     
1019                        if (lp==NULL){cout << "LP is NULL" << endl;}
1020#ifdef gfan_DEBUG
1021//                      dd_WriteLP(stdout,lp);
1022#endif 
1023                                       
1024                        lpInt=dd_MakeLPforInteriorFinding(lp);
1025                        if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
1026#ifdef gfan_DEBUG
1027//                      dd_WriteLP(stdout,lpInt);
1028#endif                 
1029
1030                        dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
1031                        if (err!=dd_NoError)
1032                        {
1033                                cout << "Error during dd_FindRelativeInterior in gcone::interiorPoint" << endl;
1034                                dd_WriteErrorMessages(stdout, err);
1035                        }
1036                       
1037                        //dd_LPSolve(lpInt,solver,&err);        //This will not result in a point from the relative interior
1038                        if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
1039                        //cout << "Tick 5" << endl;
1040                                                                       
1041                        //lpSol=dd_CopyLPSolution(lpInt);
1042                        if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}                 
1043                        //cout << "Tick 6" << endl;
1044#ifdef gfan_DEBUG
1045                        cout << "Interior point: ";
1046#endif
1047                        for (int ii=1; ii<(lpSol->d)-1;ii++)
1048                        {
1049#ifdef gfan_DEBUG
1050                                dd_WriteNumber(stdout,lpSol->sol[ii]);
1051#endif
1052                                /* NOTE This works only as long as gmp returns fractions with the same denominator*/
1053                                (iv)[ii-1]=(int)mpz_get_d(mpq_numref(lpSol->sol[ii]));  //looks evil, but does the trick
1054                        }
1055                        dd_FreeLPSolution(lpSol);
1056                        dd_FreeLPData(lpInt);
1057                        dd_FreeLPData(lp);
1058                        set_free(ddlinset);
1059                        set_free(ddredrows);                   
1060                       
1061                }//void interiorPoint(dd_MatrixPtr const &M)
1062               
1063                /** \brief Copy a ring and add a weighted ordering in first place
1064                * Kudos to walkSupport.cc
1065                */
1066                ring rCopyAndAddWeight(ring const &r, intvec const *ivw)                               
1067                {
1068                        ring res=(ring)omAllocBin(ip_sring_bin);
1069                        memcpy4(res,r,sizeof(ip_sring));
1070                        res->VarOffset = NULL;
1071                        res->ref=0;
1072                       
1073                        if (r->algring!=NULL)
1074                                r->algring->ref++;
1075                        if (r->parameter!=NULL)
1076                        {
1077                                res->minpoly=nCopy(r->minpoly);
1078                                int l=rPar(r);
1079                                res->parameter=(char **)omAlloc(l*sizeof(char_ptr));
1080                                int i;
1081                                for(i=0;i<rPar(r);i++)
1082                                {
1083                                        res->parameter[i]=omStrDup(r->parameter[i]);
1084                                }
1085                        }
1086                       
1087                        int i=rBlocks(r);
1088                        int jj;
1089                       
1090                        res->order =(int *)omAlloc((i+1)*sizeof(int));
1091                        res->block0=(int *)omAlloc((i+1)*sizeof(int));
1092                        res->block1=(int *)omAlloc((i+1)*sizeof(int));
1093                        res->wvhdl =(int **)omAlloc((i+1)*sizeof(int**));
1094                        for(jj=0;jj<i;jj++)
1095                        {                               
1096                                if (r->wvhdl[jj] != NULL)
1097                                {
1098                                        res->wvhdl[jj] = (int*) omMemDup(r->wvhdl[jj-1]);
1099                                }
1100                                else
1101                                {
1102                                        res->wvhdl[jj+1]=NULL;
1103                                }
1104                        }
1105                       
1106                        for (jj=0;jj<i;jj++)
1107                        {
1108                                res->order[jj+1]=r->order[jj];
1109                                res->block0[jj+1]=r->block0[jj];
1110                                res->block1[jj+1]=r->block1[jj];
1111                        }                                               
1112                       
1113                        res->order[0]=ringorder_a;
1114                        res->order[1]=ringorder_dp;     //basically useless, since that should never be used                   
1115                        int length=ivw->length();
1116                        int *A=(int *)omAlloc(length*sizeof(int));
1117                        for (jj=0;jj<length;jj++)
1118                        {                               
1119                                A[jj]=(*ivw)[jj];                               
1120                        }                       
1121                        res->wvhdl[0]=(int *)A;
1122                        res->block0[0]=1;
1123                        res->block1[0]=length;
1124                       
1125                        res->names = (char **)omAlloc0(rVar(r) * sizeof(char_ptr));
1126                        for (i=rVar(res)-1;i>=0; i--)
1127                        {
1128                                res->names[i] = omStrDup(r->names[i]);
1129                        }                       
1130                        rComplete(res);
1131                        return res;
1132                }//rCopyAndAdd
1133               
1134                ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
1135                {
1136                        ring res=rCopy0(currRing);
1137                        rComplete(res);
1138                        rSetWeightVec(res,(int64*)ivw);
1139                        //rChangeCurrRing(rnew);
1140                        return res;
1141                }
1142               
1143                /**
1144                * Determines whether a given facet of a cone is the search facet of a neighbouring cone
1145                * This is done in the following way:
1146                * We loop through all facets of the cone and find the "smallest" facet, i.e. the unique facet
1147                * that is first crossed during the generic walk.
1148                * We then check whether the fNormal of this facet is parallel to the fNormal of our testfacet.
1149                * If this is the case, then our facet is indeed a search facet and TRUE is retuned.
1150                */
1151                bool isSearchFacet(gcone &gcTmp, facet *testfacet)
1152                {                               
1153                        ring actRing=currRing;
1154                        facet *facetPtr=(facet*)gcTmp.facetPtr;                 
1155                        facet *fMin=new facet(*facetPtr);       //Pointer to the "minimal" facet
1156                        //facet *fMin = new facet(tmpcone.facetPtr);
1157                        //fMin=tmpcone.facetPtr;                //Initialise to first facet of tmpcone
1158                        facet *fAct;    //Ptr to alpha_i
1159                        facet *fCmp;    //Ptr to alpha_j
1160                        fAct = fMin;
1161                        fCmp = fMin->next;                             
1162                       
1163                        rChangeCurrRing(this->rootRing);        //because we compare the monomials in the rootring                     
1164                        poly p=pInit();
1165                        poly q=pInit();
1166                        intvec *alpha_i = new intvec(this->numVars);                   
1167                        intvec *alpha_j = new intvec(this->numVars);
1168                        intvec *sigma = new intvec(this->numVars);
1169                        sigma=gcTmp.getIntPoint();
1170                       
1171                        int *u=(int *)omAlloc((this->numVars+1)*sizeof(int));
1172                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
1173                        u[0]=0; v[0]=0;
1174                        int weight1,weight2;
1175                        while(fAct->next->next!=NULL)   //NOTE this is ugly. Can it be done without fCmp?
1176                        {
1177                                /* Get alpha_i and alpha_{i+1} */
1178                                alpha_i=fAct->getFacetNormal();
1179                                alpha_j=fCmp->getFacetNormal();
1180#ifdef gfan_DEBUG
1181                                alpha_i->show(); 
1182                                alpha_j->show();
1183#endif
1184                                /*Compute the dot product of sigma and alpha_{i,j}*/
1185                                weight1=dotProduct(sigma,alpha_i);
1186                                weight2=dotProduct(sigma,alpha_j);
1187#ifdef gfan_DEBUG
1188                                cout << "weight1=" << weight1 << " " << "weight2=" << weight2 << endl;
1189#endif
1190                                /*Adjust alpha_i and alpha_i+1 accordingly*/
1191                                for(int ii=1;ii<=this->numVars;ii++)
1192                                {                                       
1193                                        u[ii]=weight1*(*alpha_i)[ii-1];
1194                                        v[ii]=weight2*(*alpha_j)[ii-1];
1195                                }                               
1196                               
1197                                /*Now p_weight and q_weight need to be compared as exponent vectors*/
1198                                pSetCoeff0(p,nInit(1));
1199                                pSetCoeff0(q,nInit(1));
1200                                pSetExpV(p,u); 
1201                                pSetm(p);                       
1202                                pSetExpV(q,v); 
1203                                pSetm(q);
1204#ifdef gfan_DEBUG                               
1205                                pWrite(p);pWrite(q);
1206#endif 
1207                                /*We want to check whether x^p < x^q
1208                                => want to check for return value 1 */
1209                                if (pLmCmp(p,q)==1)     //i.e. x^q is smaller
1210                                {
1211                                        fMin=fCmp;
1212                                        fAct=fMin;
1213                                        fCmp=fCmp->next;
1214                                }
1215                                else
1216                                {
1217                                        //fAct=fAct->next;
1218                                        if(fCmp->next!=NULL)
1219                                        {
1220                                                fCmp=fCmp->next;
1221                                        }
1222                                        else
1223                                        {
1224                                                fAct=fAct->next;
1225                                        }
1226                                }
1227                                //fAct=fAct->next;
1228                        }//while(fAct.facetPtr->next!=NULL)
1229                        delete alpha_i,alpha_j,sigma;
1230                       
1231                        /*If testfacet was minimal then fMin should still point there */
1232                       
1233                        //if(fMin->getFacetNormal()==ivNeg(testfacet.getFacetNormal()))                 
1234#ifdef gfan_DEBUG
1235                        cout << "Checking for parallelity" << endl <<" fMin is";
1236                        fMin->printNormal();
1237                        cout << "testfacet is ";
1238                        testfacet->printNormal();
1239                        cout << endl;
1240#endif
1241                        if (fMin==gcTmp.facetPtr)                       
1242                        //if(areEqual(fMin->getFacetNormal(),ivNeg(testfacet.getFacetNormal())))
1243                        //if (isParallel(fMin->getFacetNormal(),testfacet->getFacetNormal()))
1244                        {                               
1245                                cout << "Parallel" << endl;
1246                                rChangeCurrRing(actRing);
1247                                //delete alpha_min, test;
1248                                return TRUE;
1249                        }
1250                        else 
1251                        {
1252                                cout << "Not parallel" << endl;
1253                                rChangeCurrRing(actRing);
1254                                //delete alpha_min, test;
1255                                return FALSE;
1256                        }
1257                }//bool isSearchFacet
1258               
1259                /** \brief Check for equality of two intvecs
1260                */
1261                bool areEqual(intvec const &a, intvec const &b)
1262                {
1263                        bool res=TRUE;
1264                        for(int ii=0;ii<this->numVars;ii++)
1265                        {
1266                                if(a[ii]!=b[ii])
1267                                {
1268                                        res=FALSE;
1269                                        break;
1270                                }
1271                        }
1272                        return res;
1273                }
1274               
1275                /** \brief The reverse search algorithm
1276                */
1277                void reverseSearch(gcone *gcAct) //no const possible here since we call gcAct->flip
1278                {
1279                        facet *fAct=new facet();
1280                        fAct = gcAct->facetPtr;                 
1281                       
1282                        while(fAct->next!=NULL)  //NOTE NOT SURE WHETHER THIS IS RIGHT! Do I reach EVERY facet or only all but the last?
1283                        {
1284                                cout << "==========================================================================================="<< endl;
1285                                gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
1286                                gcone *gcTmp = new gcone(*gcAct);
1287                                //idShow(gcTmp->gcBasis);
1288                                gcTmp->getConeNormals(gcTmp->gcBasis, TRUE);
1289#ifdef gfan_DEBUG
1290                                facet *f = new facet();
1291                                f=gcTmp->facetPtr;
1292                                while(f->next!=NULL)
1293                                {
1294                                        f->printNormal();
1295                                        f=f->next;                                     
1296                                }
1297#endif
1298                                gcTmp->showIntPoint();
1299                                /*recursive part goes gere*/
1300                                if (isSearchFacet(*gcTmp,(facet*)gcAct->facetPtr))
1301                                {
1302                                        gcAct->next=gcTmp;
1303                                        cout << "PING"<< endl;
1304                                        reverseSearch(gcTmp);
1305                                }
1306                                else
1307                                {
1308                                        delete gcTmp;
1309                                        /*NOTE remove fAct from linked list. It's no longer needed*/
1310                                }
1311                                /*recursion ends*/
1312                                fAct = fAct->next;             
1313                        }//while(fAct->next!=NULL)
1314                }//reverseSearch
1315               
1316                /** \brief The new method of Markwig and Jensen
1317                */
1318                void noRevS(gcone &gc)
1319                {
1320                        facet *fListPtr = new facet();                 
1321                        facet *fAct = new facet();
1322                        fAct = gc.facetPtr;
1323                                               
1324                        dd_set_global_constants();
1325                        dd_rowrange ddrows;
1326                        dd_colrange ddcols;
1327                        ddrows=2;       //Each facet is described by its normal
1328                        ddcols=gc.numVars+1;    // plus one for the standard simplex
1329                       
1330                        while(fAct->next!=NULL)
1331                        {
1332                                /*Compute an interior point for each facet*/                           
1333                                dd_MatrixPtr ddineq;   
1334                                ddineq=dd_CreateMatrix(ddrows,ddcols);
1335                                intvec *comp = new intvec(this->numVars);
1336                                comp=fAct->getFacetNormal();                           
1337                                for(int ii=0; ii<this->numVars; ii++)                                   
1338                                {                                                                               
1339                                        dd_set_si(ddineq->matrix[0][ii+1],(*comp)[ii]);
1340                                        dd_set_si(ddineq->matrix[1][ii+1],1);   //Assure we search in the pos. orthant                 
1341                                }
1342                                set_addelem(ddineq->linset,1);  //We want equality in the first row
1343                                //dd_WriteMatrix(stdout,ddineq);
1344                                interiorPoint(ddineq,*comp);                           
1345                                /**/
1346#ifdef gfan_DEBUG
1347                                cout << "IP is";
1348                                comp->show(); cout << endl;
1349#endif
1350                                //Store the interior point and the UCN
1351                                fListPtr->setInteriorPoint( comp );                             
1352                                fListPtr->setUCN( gc.getUCN() );       
1353                                                       
1354                                dd_FreeMatrix(ddineq);
1355                                fAct=fAct->next;        //iterate
1356                        }                       
1357                       
1358                        //NOTE Hm, comment in and get a crash for free...
1359                        //dd_free_global_constants();                           
1360                        gc.writeConeToFile(gc);
1361                       
1362                        /*2nd step
1363                        Choose a facet from fListPtr, flip it and forget the previous cone
1364                        */
1365                        fAct = fListPtr;
1366                        gcone *gcTmp = new gcone(gc);   //copy
1367                        gcTmp->flip(gcTmp->gcBasis,fAct);
1368                        //NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
1369                       
1370                }//void noRevS(gcone &gc)
1371               
1372                /** \brief Write information about a cone into a file on disk
1373                *
1374                * This methods writes the information needed for the "second" method into a file.
1375                * The file's is divided in sections containing information on
1376                * 1) the ring
1377                * 2) the cone's Groebner Basis
1378                * 3) the cone's facets
1379                * Each line contains exactly one date
1380                * Each section starts with its name in CAPITALS
1381                */
1382                void writeConeToFile(gcone const &gc)
1383                {
1384                        ofstream gcOutputFile("/tmp/cone1.gc");
1385                        if (!gcOutputFile)
1386                        {
1387                                cout << "Error opening file for writing in writeConeToFile" << endl;
1388                        }
1389                        else
1390                        {                       
1391                                gcOutputFile << "RING" << endl;                         
1392                                //Write this->gcBasis into file
1393                                gcOutputFile << "GCBASIS" << endl;                             
1394                                for (int ii=0;ii<IDELEMS(gc.gcBasis);ii++)
1395                                {                                       
1396                                        gcOutputFile << p_String((poly)gc.gcBasis->m[ii],gc.baseRing) << endl;
1397                                }                       
1398                               
1399                                facet *fAct = new facet();
1400                                fAct = gc.facetPtr;
1401                                gcOutputFile << "FACETS" << endl;                                                               
1402                                while(fAct!=NULL)
1403                                {       
1404                                        intvec *iv = new intvec(gc.numVars);
1405                                        iv=fAct->getFacetNormal();
1406                                        for (int ii=0;ii<iv->length();ii++)
1407                                        {
1408                                                if (ii<iv->length()-1)
1409                                                {
1410                                                        gcOutputFile << (*iv)[ii] << ",";
1411                                                }
1412                                                else
1413                                                {
1414                                                        gcOutputFile << (*iv)[ii] << endl;
1415                                                }
1416                                        }
1417                                        fAct=fAct->next;
1418                                }                               
1419                               
1420                                gcOutputFile.close();
1421                                delete fAct; fAct=NULL;
1422                        }
1423                }//writeConeToFile(gcone const &gc)
1424               
1425                /** \brief Reads a cone from a file identified by its number
1426                */
1427                void readConeFromFile(int gcNum)
1428                {
1429                }
1430               
1431        friend class facet;
1432};//class gcone
1433
1434ideal gfan(ideal inputIdeal)
1435{
1436        int numvar = pVariables; 
1437       
1438        enum searchMethod {
1439                reverseSearch,
1440                noRevS
1441        };
1442       
1443        searchMethod method;
1444        method = noRevS;
1445//      method = reverseSearch;
1446       
1447#ifdef gfan_DEBUG
1448        cout << "Now in subroutine gfan" << endl;
1449#endif
1450        ring inputRing=currRing;        // The ring the user entered
1451        ring rootRing;                  // The ring associated to the target ordering
1452        ideal res;     
1453        facet *fRoot;
1454       
1455        if (method==reverseSearch)
1456        {
1457       
1458        /* Construct a new ring which will serve as our root*/
1459        rootRing=rCopy0(currRing);
1460        rootRing->order[0]=ringorder_lp;
1461       
1462        rComplete(rootRing);
1463        rChangeCurrRing(rootRing);
1464       
1465        /* Fetch the inputIdeal into our rootRing */
1466        map theMap=(map)idMaxIdeal(1);  //evil hack!
1467        theMap->preimage=NULL;  //neccessary?
1468        ideal rootIdeal;
1469        rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
1470#ifdef gfan_DEBUG
1471        cout << "Root ideal is " << endl;
1472        idShow(rootIdeal);
1473        cout << "The root ring is " << endl;
1474        rWrite(rootRing);
1475        cout << endl;
1476#endif 
1477       
1478        //gcone *gcRoot = new gcone();  //Instantiate the sink
1479        gcone *gcRoot = new gcone(rootRing,rootIdeal);
1480        gcone *gcAct;
1481        gcAct = gcRoot;
1482        gcAct->numVars=pVariables;
1483        gcAct->getGB(rootIdeal);        //sets gcone::gcBasis
1484        idShow(gcAct->gcBasis);
1485        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
1486        //gcAct->flip(gcAct->gcBasis,gcAct->facetPtr); 
1487        /*Now it is time to compute the search facets, respectively start the reverse search.
1488        But since we are in the root all facets should be search facets. IS THIS TRUE?
1489        NOTE: Check for flippability is not very sophisticated
1490        */     
1491        //gcAct->reverseSearch(gcAct); 
1492        rChangeCurrRing(rootRing);
1493        res=gcRoot->gcBasis;   
1494        }//if method==reverSearch
1495       
1496        if(method==noRevS)
1497        {
1498                gcone *gcRoot = new gcone(currRing,inputIdeal);
1499                gcone *gcAct;
1500                gcAct = gcRoot;
1501                gcAct->numVars=pVariables;
1502                gcAct->getGB(inputIdeal);
1503                gcAct->getConeNormals(gcAct->gcBasis);
1504                cout << "ding" << endl;         
1505                gcAct->noRevS(*gcAct);         
1506                res=gcAct->gcBasis;     
1507        }
1508       
1509        /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
1510        The return type will then be of type LIST_CMD
1511        Assume gfan has finished, thus we have enumerated all the cones
1512        Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
1513        Groebner Basis and merge this somehow with LIST_CMD
1514        => Count the cones!
1515        */
1516        //rChangeCurrRing(rootRing);
1517        //res=gcAct->gcBasis;
1518        //res=gcRoot->gcBasis; 
1519        return res;
1520        //return GBlist;
1521}
1522/*
1523Since gfan.cc is #included from extra.cc there must not be a int main(){} here
1524*/
1525#endif
Note: See TracBrowser for help on using the repository browser.