Changeset 5cd07b in git


Ignore:
Timestamp:
Apr 3, 2009, 3:49:16 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')
Children:
509339932269ee5ca46faa0b1aa9cfd1de39df76
Parents:
c4a041219a2ad7d07a655aa5c8e89a483611a985
Message:
Two new methods:
int gcone::dotProduct(intvec **a, intvec **b)
bool gcone::isParallel(intvec *a, intvec *b)
gcone::flip now computes the initial form as an array of polynomials.
How to get this into an ideal?


git-svn-id: file:///usr/local/Singular/svn/trunk@11615 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
kernel
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    rc4a041 r5cd07b  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-04-02 09:44:23 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.27 2009-04-02 09:44:23 monerjan Exp $
    6 $Id: gfan.cc,v 1.27 2009-04-02 09:44:23 monerjan Exp $
     4$Date: 2009-04-03 13:49:16 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.28 2009-04-03 13:49:16 monerjan Exp $
     6$Id: gfan.cc,v 1.28 2009-04-03 13:49:16 monerjan Exp $
    77*/
    88
    … …  
    281281                                /*Now load should be full and we can call setFacetNormal*/
    282282                                fAct->setFacetNormal(load);
    283                                 fAct->printNormal();
     283                                //fAct->printNormal();
    284284                                fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
    285285                                delete load;
    … …  
    301301                }//method getConeNormals(ideal I)       
    302302               
    303                 bool isParallel(int v[], intvec iv)
    304                 {
    305                 }
     303                /*bool isParallel(int v[], intvec iv)
     304                {
     305}*/
    306306               
    307307                /** \brief Compute the Groebner Basis on the other side of a shared facet
    … …  
    311311                * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
    312312                * is parallel to \f$ leadexp(g) \f$
    313                 * Checking for parallelity is done by brute force dividing of components.
    314                 * Other possibilitiesincludes  computing the rank of the matrix consisting of the vectors in question and
     313                * Parallelity is checked using basic linear algebra. See gcone::isParallel.
     314                * Other possibilities includes  computing the rank of the matrix consisting of the vectors in question and
    315315                * computing an interior point of the facet and taking all terms having the same weight with respect
    316316                * to this interior point.
    … …  
    322322                       
    323323                        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
    324                         fNormal = f->getFacetNormal();  //read this->fNormal;                   
     324                        fNormal = f->getFacetNormal();  //read this->fNormal;
     325#ifdef gfan_DEBUG
     326                        cout << "fNormal=";
    325327                        fNormal->show();
    326                        
    327                         bool isParallel=TRUE;
    328                         double lambda=0;
    329                         double lambda_temp=lambda;
    330                         /* EXTREMELY EXPERIMENTAL CODE AHEAD - NUMERICAL INSTABILITIES GUARANTEED DUE TO DIVISONS*/
     328                        cout << endl;
     329#endif                         
    331330                        /*1st step: Compute the initial ideal*/
    332                         poly initialForm[IDELEMS(gb)];  //array of #polys in GB to store initial form
     331                        poly initialFormElement[IDELEMS(gb)];   //array of #polys in GB to store initial form
     332                        ideal initialForm;
    333333                        poly aktpoly;
    334                         int check[this->numVars];       //array to store the difference of LE and v
     334                        intvec *check = new intvec(this->numVars);      //array to store the difference of LE and v
    335335                       
    336336                        for (int ii=0;ii<IDELEMS(gb);ii++)
    … …  
    338338                                aktpoly = (poly)gb->m[ii];                                                             
    339339                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    340                                 int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));                           
    341                                 intvec *iv=(intvec *)omAlloc((this->numVars+1)*sizeof(intvec)); //let's try this
    342                                 intvec *ivLeadExpV=(intvec *)omAlloc((this->numVars+1)*sizeof(intvec));
     340                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
    343341                                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
    344342                                //pGetExpV(aktpoly,ivLeadExpV);
    345                                 initialForm[ii]=pHead(aktpoly);
     343                                initialFormElement[ii]=pHead(aktpoly);
    346344                               
    347345                                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
    348346                                {
    349347                                        aktpoly=pNext(aktpoly); //next term
    350                                         pGetExpV(aktpoly,v);
    351                                        
    352                                         for (int i=0;i<pVariables;i++)
     348                                        pSetm(aktpoly);
     349                                        pGetExpV(aktpoly,v);           
     350                                        /* Convert (int)v into (intvec)check */                 
     351                                        for (int jj=0;jj<this->numVars;jj++)
    353352                                        {
    354                                                 check[i] = v[i]-leadExpV[i];
    355                                                 cout << "check["<<i<<"]="<<check[i]<<endl;
     353                                                //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
     354                                                //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
     355                                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
    356356                                        }
    357                                         lambda=check[0]/(*fNormal)[0];
    358                                         lambda_temp=lambda;
    359                                         cout << "lambda="<<lambda<<endl;
    360                                         for (int i=0;i<pVariables;i++)
     357                                        cout << "check=";                       
     358                                        check->show();
     359                                        cout << endl;
     360                                        if (isParallel(check,fNormal)) //pass *check when
    361361                                        {
    362                                                 lambda_temp=check[i]/(*fNormal)[i];
    363                                                 cout<<"lambda_temp="<<lambda_temp<<endl;
    364                                                 if (lambda_temp!=lambda)
    365                                                 {
    366                                                         isParallel=FALSE;
    367                                                 }
    368                                         }
    369                                         if (isParallel==TRUE) //pass *check when
    370                                         {
    371                                                 cout << "PARALLEL" << endl;
    372                                                 //initialForm[ii] += pHead(aktpoly);    //This should add the respective term to
    373                                         }                               
    374                                 }
    375                                
    376                         }
     362                                                cout << "Parallel vector found, adding to initialFormElement" << endl;                         
     363                                                initialFormElement[ii] = pAdd(initialFormElement[ii],(poly)pHead(aktpoly));
     364                                        }                                               
     365                                }//while
     366                                cout << "Initial Form=";                               
     367                                pWrite(initialFormElement[ii]);
     368                                cout << "---" << endl;
     369                                /*Now initialFormElement must be added to (ideal)initialForm */                                                 
     370                        }//for
    377371                }//void flip(ideal gb, facet *f)
    378372                               
    … …  
    392386               
    393387                ideal GenGrbWlk(ideal, ideal);  //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
    394                
    395 
     388
     389
     390                /** \brief Compute the dot product of two intvecs
     391                *
     392                * THIS IS WEIRD - BUT WORKS
     393                */
     394                int dotProduct(intvec **a, intvec **b)                         
     395                {
     396                        intvec iva=*a;
     397                        intvec ivb=*b;
     398                        int res=0;
     399                        for (int i=0;i<this->numVars;i++)
     400                        {
     401                                res = res+(iva[i]*ivb[i]);
     402                        }
     403                        return res;
     404                }
     405
     406                /** \brief Check whether two intvecs are parallel
     407                *
     408                * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
     409                */
     410                bool isParallel(intvec *a, intvec *b)
     411                {                       
     412                        //bool res=FALSE;
     413                        //intvec iva(this->numVars);
     414                        //intvec ivb(this->numVars);
     415                        int lhs,rhs;
     416                        //lhs=0;
     417                        //rhs=0;
     418                        //iva = a;
     419                        //ivb = b;
     420                        //lhs=dotProduct(&iva,&ivb)*dotProduct(&iva,&ivb);
     421                        //lhs=dotProduct(&iva,&iva)*dotProduct(&ivb,&ivb);
     422                        lhs=dotProduct(&a,&b)*dotProduct(&a,&b);
     423                        rhs=dotProduct(&a,&a)*dotProduct(&b,&b);
     424                        cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
     425                        if (lhs==rhs)
     426                        {
     427                                return TRUE;
     428                        }
     429                        else
     430                        {
     431                                return FALSE;
     432                        }
     433                }//bool isParallel
     434               
     435                /** \brief Check two intvecs for equality --- PROBABLY NOT NEEDED
     436                *
     437                * \f$ \alpha=\beta\Leftrightarrow\langle\alpha-\beta,\alpha-\beta\rangle=0 \f$
     438                */
     439                bool isEqual(intvec a, intvec b)
     440                {
     441                        intvec *ivdiff;
     442                        int res;
     443                       
     444                        ivdiff=ivSub(&a,&b);
     445                        cout << "ivdiff=";
     446                        ivdiff->show();
     447                        cout << endl;
     448                        //res=dotProduct(ivdiff,ivdiff);
     449                        if (res==0)
     450                        {
     451                                return TRUE;
     452                        }
     453                        else
     454                        {
     455                                return FALSE;
     456                        }                       
     457                }//bool isEqual
    396458};//class gcone
    397 
    398 /*
    399 function getGB incorporated into class gcone with rev 1.24
    400 */
    401 
    402 //DEPRECATED since rev 1.24 with existence of gcone::getConeNormals(ideal I);
    403 //Kept for unknown reasons ;)
    404 facet *getConeNormals(ideal I)
    405 {
    406         return NULL;
    407 }
    408459
    409460ideal gfan(ideal inputIdeal)
    … …  
    433484        rootRing->order[0]=ringorder_dp;
    434485        rComplete(rootRing);
    435         rChangeCurrRing(rootRing);
     486        //rChangeCurrRing(rootRing);
    436487        ideal rootIdeal;
    437488        /* Fetch the inputIdeal into our rootRing */
    438489        map m=(map)idInit(IDELEMS(inputIdeal),0);
    439         rootIdeal=fast_map(inputIdeal,inputRing,(ideal)m, currRing);
     490        //rootIdeal=fast_map(inputIdeal,inputRing,(ideal)m, currRing);
    440491#ifdef gfan_DEBUG
    441492        cout << "Root ideal is " << endl;

  • kernel/gfan.h

    rc4a041 r5cd07b  
    22gfan.h Interface to gfan.cc
    33
    4 $Author: Singular $
    5 $Date: 2009-02-19 15:39:08 $
    6 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.h,v 1.5 2009-02-19 15:39:08 Singular Exp $
    7 $Id: gfan.h,v 1.5 2009-02-19 15:39:08 Singular Exp $
     4$Author: monerjan $
     5$Date: 2009-04-03 13:49:16 $
     6$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.h,v 1.6 2009-04-03 13:49:16 monerjan Exp $
     7$Id: gfan.h,v 1.6 2009-04-03 13:49:16 monerjan Exp $
    88*/
    99#ifdef HAVE_GFAN
    1010#include "intvec.h"
    1111
    12 ideal getGB(ideal inputIdeal);
     12//ideal getGB(ideal inputIdeal);
    1313ideal gfan(ideal inputIdeal);
     14//int dotProduct(intvec a, intvec b);
     15//bool isParallel(intvec a, intvec b);
    1416#endif
Note: See TracChangeset for help on using the changeset viewer.