=================================================================== RCS file: /home/cvs/OpenXM/src/kan96xx/Doc/gfan.sm1,v retrieving revision 1.9 retrieving revision 1.17 diff -u -p -r1.9 -r1.17 --- OpenXM/src/kan96xx/Doc/gfan.sm1 2005/06/30 08:39:39 1.9 +++ OpenXM/src/kan96xx/Doc/gfan.sm1 2009/09/04 11:13:11 1.17 @@ -1,9 +1,25 @@ -% $OpenXM: OpenXM/src/kan96xx/Doc/gfan.sm1,v 1.8 2004/10/13 23:36:52 takayama Exp $ +% $OpenXM: OpenXM/src/kan96xx/Doc/gfan.sm1,v 1.16 2009/09/04 02:59:55 takayama Exp $ % cp cone.sm1 $OpenXM_HOME/src/kan96xx/Doc/gfan.sm1 -% $Id: gfan.sm1,v 1.9 2005/06/30 08:39:39 takayama Exp $ +% $Id: gfan.sm1,v 1.17 2009/09/04 11:13:11 takayama Exp $ % iso-2022-jp %%Ref: @s/2004/08/21-note.pdf +%% gfan.sm1 works only for polymake 2.0 Use webservice of 2.0. +[(gfan) +[ + (gfan.sm1 is a package to compute global and local Grobner fans.) + (See R.Bahloul and N.Takayama, arxiv, math.AG/0412044 and references as to algorithms.) + (At the beginning of gfan.sm1, there are sample inputs cone.sample and cone.sample2.) + ( ) + (gfan.sm1 works only with polymake 2.0. We provides a web service of computing ) + (with polymake 2.0. /@@@polymake.web 1 def is set by default in gfan.sm1.) + (See changelog-ja.tex as to details on the difference between 2.0 and later versions.) + ( ) + ( cone.Wt cone.Lpt {vertices in the output} are weights on the rays of the Grobner cone.) + ( cone.L gives a basis of the linearity space.) +] +] putUsages + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Two examples are given below to get a global Grobner fan and %% a local Grobner fan ; cone.sample and cone.sample2 @@ -17,6 +33,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /cone.sample { cone.load.cohom + /cone.ckmFlip 1 def % write a comment about the problem. "nl" means new line. /cone.comment [ (Toric ideal for 1-simplex x 2-simplex, in k[x]) nl @@ -76,6 +93,7 @@ def ] def +/cone.DhH 0 def % Set a function to compute Grobner basis. % cone.gb_Dh : For computing in Homogenized Weyl algebra h[1,1](D). % cone.gb_DhH : For computing in doubly homogenized Weyl algebra. @@ -128,6 +146,7 @@ printGrobnerFan %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /cone.sample2 { cone.load.cohom + /cone.ckmFlip 1 def % write a comment about the problem. "nl" means new line. /cone.comment [ (BS for y and y-(x-1)^2, t1, t2 space, in doubly homogenized Weyl algebra.) nl @@ -194,6 +213,7 @@ def cone.input { . homogenize toString } map /cone.input set dh.end +/cone.DhH 1 def % Set a function to compute Grobner basis. % cone.gb_Dh : For computing in Homogenized Weyl algebra h[1,1](D). % cone.gb_DhH : For computing in doubly homogenized Weyl algebra. @@ -246,9 +266,13 @@ dhcone.printGrobnerFan % If you use the cgi/polymake on the net, then uncomment out. %/doPolymake {doPolymake.OoHG} def (Using doPolymake.OoHG ) message %/polymake.start {polymake.start.OoHG} def (Using polymake.start.OoHG ) message +/@@@polymake.web 1 def %% Choose it automatically. -[(which) (polymake)] oxshell tag 0 eq { - (Polymake is not installed in this system.) message +[(which) (polymake)] oxshell tag 0 eq +@@@polymake.web 1 eq +or +{ + (Polymake is not installed in this system or @@@polymake.web is set.) message /doPolymake {doPolymake.OoHG} def (Using doPolymake.OoHG ) message /polymake.start {polymake.start.OoHG} def @@ -398,8 +422,36 @@ cone.comment message % cone.d (pointed cones lies in this space. cf. cone.Lp) % These are set during getting the cone.startingCone +%< +% global +%cone.ckmFlip. Collar-Kalkbrener-Mall の flip アルゴリズムを使わない 0. 使う 1. +% Default は 0. +%> +/cone.ckmFlip 0 def %< +% global +% cone.DhH dx x = x dx + h H なら 1. dx x = x dx + h^2 なら 0. Default 0. +%> +/cone.DhH 0 def + +%< +% Global +% gbCheck をするか? しないと結果はあやふや. しかしメモリ exhaust は防げる. +% 使うときは /cone.epsilon, /cone.epsilon.limit を十分小さくしておく. +%> +/cone.do_gbCheck 1 def + +% Default の cone.gb の定義. 各プログラムで再度定義してもよい. +/cone.gb { + cone.DhH { + cone.gb_DhH + } { + cone.gb_Dh + } ifelse +} def + +%< % Usage: wv g coneEq1 % in(f) が monomial 専用. in_w(f) = LT(f) となる weight w の満たす % 不等式制約を求める. @@ -988,17 +1040,31 @@ def %> /getConeInfo { /arg1 set - [/ww /g /ceq /ceq2 /cdim /mmc /mmL /rr /ineq /ppt] pushVariables + [/ww /g /ceq /ceq2 /cdim /mmc /mmL /rr /ineq /ppt /rr0 /mm0 /mm1] pushVariables [ /ceq arg1 def ceq pruneZeroVector /ceq set + + ceq length 0 eq { + (Monomial ideal is not accepted as an input.) cone_ir_input + } { } ifelse + + /mm1 + ( Use [(keep_tmp_files) 1] oxshell to check the input to polymake2tfb. See /tmp or $TMP ) + def + ceq genPo2 /ceq2 set % ceq2 は polymake.data(polymake.INEQUALITIES(...)) 形式 % polymake で ceq2 の次元の計算. /getConeInfo.ceq ceq def /getConeInfo.ceq2 ceq2 def cone.debug { (Calling polymake DIM.) message } { } ifelse - [(DIM) ceq2] doPolymake 1 get /rr set + [(DIM) ceq2] doPolymake /rr0 set + % rr0 2 get message + rr0 2 get 1 get 0 get /mm0 set + mm0 length 0 eq { } + { [mm0 mm1] cat error } ifelse + rr0 1 get /rr set cone.debug {(Done.) message } { } ifelse % test5 には次のコメントとりさる. 上の行をコメントアウト. % test5.data tfbToTree /rr set @@ -1013,13 +1079,24 @@ def % FACETS を持っていないなら再度計算する. % POINTED, NOT__POINTED も得られる cone.debug { (Calling polymake FACETS.) message } { } ifelse - [(FACETS) ceq2] doPolymake 1 get /rr set + [(FACETS) ceq2] doPolymake /rr0 set + + % rr0 2 get message + rr0 2 get 1 get 0 get /mm0 set + mm0 length 0 eq { } + { [mm0 mm1] cat error } ifelse + + rr0 1 get /rr set cone.debug { (Done.) message } { } ifelse } { } ifelse rr (VERTICES) getNode tag 0 eq { (internal error: VERTICES is not found.) error - } { } ifelse + } { + rr (VERTICES) getNode + (UNDEF) getNode tag 0 eq { } + { (internal error: VERTICES is UNDEF. See rr. Set /@@@polymake.web 1 def) error } ifelse + } ifelse /cone.getConeInfo.rr1 rr def @@ -1589,7 +1666,7 @@ def /vlist arg1 def /wlist arg2 def wlist length vlist length eq { - } { (cone_wtowv: length of the argument must be the same.) error} ifelse + } { (cone_wtowv: length of the argument must be the same. Please check the values of cone.vlist cone.vv cone.type parametrizeWeightSpace) error} ifelse wlist to_int32 /wlist set [ @@ -1846,6 +1923,11 @@ def % note: 2004.9.2 cone (facetsv) getNode 2 get facet_i get /v set cone (facets) getNode 2 get facet_i get /f set + + v length 0 eq { + (The codimension of the linarity space of the Grobner cone seems to be 1 or 0.) cone_ir_input + } { } ifelse + /vp v 0 get def 1 1 v length 1 sub { /ii set @@ -2023,15 +2105,20 @@ def %> /cone.gb_Dh { /arg2 set /arg1 set - [/ff /ww /gg] pushVariables + [/ff /ww /gg /gbopt] pushVariables [ /ff arg1 def /ww arg2 def [(AutoReduce) 1] system_variable [cone.vv ring_of_differential_operators [ww] weight_vector 0] define_ring - [ff {toString .} map] ff getAttributeList setAttributeList - groebner 0 get /gg set + %(---) messagen ff getAttributeList message + ff getAttributeList tag 0 eq {/gbopt [ ] def } + { + /gbopt ff getAttributeList def + } ifelse + [ff {toString .} map gbopt] + groebner 0 get /gg set %% groenber は attribute を受け付けない. /cone.gb_Dh.g gg def /arg1 gg def ] pop @@ -2361,7 +2448,7 @@ def % Usages: result_getNextFlip getNextCone ncone % flip して新しい ncone を得る. %> -/getNextCone { +/getNextCone.orig { /arg1 set [/ncone /ccone /kk /w /next_weight_w_wv] pushVariables [ @@ -3430,6 +3517,7 @@ def % gb の check をやるので, それに失敗したら null を戻す. % weight はすべて vw 形式で. vw 形式 = variable weight の繰り返しの形式 % reducedGb は文字列のリストではなく多項式の形式のこと. +% 理由は reducedGb より ring の構造を読むため. %> /ckmFlip { /arg1 set @@ -3476,6 +3564,12 @@ def % 例. [(x) -1 (Dx) 1 (y) -1 (Dy) 2] ==> [(x) -1 (Dx) 1 (y) -1 (y') 2] facetWeight { dup isString { . rTable replace toString } { } ifelse } map /facetWeight_gr set + +% newWeight も 新しい環 gr_ww の weight に変換. +% 例. [(x) -1 (Dx) 1 (y) -1 (Dy) 2] ==> [(x) -1 (Dx) 1 (y) -1 (y') 2] + newWeight { dup isString { . rTable replace toString } + { } ifelse } map /newWeight_gr set + % Dx x = x Dx + h H or Dx x = x Dx + h^2 で計算. % どちらをとるかは cone.gb_gr で区別するしかなし %% [ch1 vlist_gr oldWeight_gr] /ttt set @@ -3484,10 +3578,6 @@ def ch1 {toString .} map /ch1 set %% ここまででとりあえずテストをしよう. %% ch1 /arg1 set -% newWeight も 新しい環 gr_ww の weight に変換. -% 例. [(x) -1 (Dx) 1 (y) -1 (Dy) 2] ==> [(x) -1 (Dx) 1 (y) -1 (y') 2] - newWeight { dup isString { . rTable replace toString } - { } ifelse } map /newWeight_gr set [ch1 { toString } map vlist_gr newWeight_gr] cone.gb_gr /ch2 set % Dx x = x Dx + h H or Dx x = x Dx + h^2 で計算. @@ -3505,15 +3595,23 @@ def ccf { 2 get {toString . rTable2 replace toString} map } map /ccf2 set %% ccf2 は gr でない ring の元. gOld getRing ring_def - cone.beginH % Hh か h^2 か. + cone.DhH { cone.begin_DhH } { } ifelse % Hh か h^2 か. ccf2 { {.} map gOld mul } map /gNew set gNew { toString } map /gNew set - cone.endH + cone.DhH { cone.end_DhH } { } ifelse % Hh か h^2 か. % gNew /arg1 set %gNew が newWeight での GB か check. Yes なら reduced basis へ. %No なら null を戻す. - gNew [(gbCheck) 1] setAttributeList newWeight - cone.gb (gb) getAttribute +%%Ref: note @s/2005/06/30-note-gfan.pdf + cone.do_gbCheck not { + (Warning! gbCheck is skipped.) message + } { + (Doing gbCheck.) message + } ifelse + cone.do_gbCheck { + gNew [(gbCheck) 1] setAttributeList newWeight + cone.gb (gb) getAttribute + } { 1 } ifelse 1 eq { gNew [(reduceOnly) 1] setAttributeList newWeight cone.gb /arg1 set }{ /arg1 null def } ifelse @@ -3524,30 +3622,122 @@ def %< % Usages: f gbasis cone.reduction_DhH +% dx x = x dx + h H での reduction. %> /cone.reduction_DhH { /arg2 set /arg1 set [/ff /ggbasis /eenv /ans] pushVariables [ /ff arg1 def /ggbasis arg2 def - cone.beginH + cone.begin_DhH ff ggbasis reduction /ans set - cone.endH + cone.end_DhH /arg1 ans def ] pop popVariables arg1 } def +%< +% Usages: f gbasis cone.reduction_Dh +% dx x = x dx + h^2 での reduction. +%> +/cone.reduction_Dh { + /arg2 set /arg1 set + [/ff /ggbasis /eenv /ans] pushVariables + [ + /ff arg1 def /ggbasis arg2 def + ff ggbasis reduction /ans set + /arg1 ans def + ] pop + popVariables + arg1 +} def + +%< +% Usages: cone.begin_DhH dx x = x dx + h H を開始. +%> /cone.begin_DhH { [(Homogenize) (AutoReduce) (KanGBmessage)] pushEnv /cone.eenv set [(Homogenize) 3] system_variable } def +%< +% Usages: cone.begin_DhH dx x = x dx + h H を終了. +%> /cone.end_DhH { cone.eenv popEnv } def +%< +% Usages: ff vv ww cone.gb_gr_DhH dx x = x dx + h H で計算. +% dh.gb は dhecart.sm1 で定義されており, dx x = x dx + h H での計算. +% gr をとっても, -w,w の場合は 微分作用素環のままであり, これが必要. +% bug? cone.gb で十分? +%> +/cone.gb_gr_DhH { + /arg1 set + [/ff /ww /vv] pushVariables + [ + /ff arg1 0 get def + /vv arg1 1 get def + /ww arg1 2 get def + /dh.gb.verbose 1 def + /dh.autoHomogenize 0 def + [(AutoReduce) 1] system_variable + [ff { toString } map vv + [ww vv generateD1_1]] dh.gb 0 get /arg1 set + ] pop + popVariables + arg1 +} def +%< +% Usages: ff vv ww cone.gb_gr_Dh dx x = x dx + h^2 で計算. +% gb は dhecart.sm1 で定義されており, dx x = x dx + h^2 での計算. +% gr をとっても, -w,w の場合は 微分作用素環のままであり, これが必要. +% bug? cone.gb で十分? +%> +/cone.gb_gr_Dh { + /arg1 set + [/ff /ww /vv /gg /envtmp] pushVariables + [ + /ff arg1 0 get def + /vv arg1 1 get def + /ww arg1 2 get def + + [(AutoReduce) (KanGBmessage)] pushEnv /envtmp set + [(AutoReduce) 1] system_variable + [(KanGBmessage) 1] system_variable + [vv ring_of_differential_operators + [ww] weight_vector 0] define_ring + [ff {toString .} map] ff getAttributeList setAttributeList + groebner 0 get /gg set + envtmp popEnv + + /arg1 gg def + ] pop + popVariables + arg1 +} def + + +% これらは cone.ckmFlip 1 の時しか使わず. +/cone.reduction { + cone.DhH { + cone.reduction_DhH + }{ + cone.reduction_Dh + } ifelse +} def +/cone.gb_gr { + cone.DhH { + cone.gb_gr_DhH + }{ + cone.gb_gr_Dh + } ifelse +} def + + /test1.ckmFlip { % cf. cone.sample2 cone.load.cohom @@ -3561,6 +3751,10 @@ def /cone.parametrizeWeightSpace { 4 2 parametrizeSmallFan } def + + /cone.DhH 1 def + /cone.ckmFlip 1 def + /cone.local 1 def /cone.w_start null def /cone.h0 1 def @@ -3579,38 +3773,8 @@ def cone.input { . homogenize toString } map /cone.input set dh.end - /cone.gb { - cone.gb_DhH - } def - /cone.reduction { - cone.reduction_DhH - } def - - /cone.beginH { - cone.begin_DhH - } def - /cone.endH { - cone.end_DhH - } def % テストを開始する. - /cone.gb_gr { - /arg1 set - [/ff /ww /vv] pushVariables - [ - /ff arg1 0 get def - /vv arg1 1 get def - /ww arg1 2 get def - /dh.gb.verbose 1 def - /dh.autoHomogenize 0 def - [(AutoReduce) 1] system_variable - [ff { toString } map vv - [ww vv generateD1_1]] dh.gb 0 get /arg1 set - ] pop - popVariables - arg1 - } def - % getStartingCone /cone.ncone set % cone.ncone updateFan % cone.gblist 0 get message @@ -3625,4 +3789,125 @@ def [ff (t1,t2,x,y) wOld wFacet wNew] ckmFlip /ff2 set (See ff and ff2) message +} def + +%< +% Usages: cone i getaVectorOnFacet +% cone の i 番目の facet の上の vector を求める. +% cf. liftWeight +%> +/getaVectorOnFacet { + /arg2 set /arg1 set + [/cone /facet_i /ep /vp /v /v /ii] pushVariables + [ + /cone arg1 def /facet_i arg2 def + facet_i to_int32 /facet_i set + + cone (facetsv) getNode 2 get facet_i get /v set + /vp v 0 get def + 1 1 v length 1 sub { + /ii set + vp v ii get add /vp set + } for + vp nnormalize_vec /vp set + /arg1 vp def + ] pop + popVariables + arg1 +} def + +/getNextCone { + getNextCone_ckm +} def + +%< +% Usages: result_getNextFlip getNextCone_ckm ncone +% flip して新しい ncone を得る. Collar-Kalkbrener-Moll のアルゴリズムを使う +% if (cone.ckmFlip == 0) 普通の計算 else CKM. +%> +/getNextCone_ckm { + /arg1 set + [/ncone /ccone /kk /w /next_weight_w_wv /cid /ttt] pushVariables + [ + /ccone arg1 def + /ncone null def + /kk ccone 1 get def % kk は cid 番目の cone の kk 番目の facet を表す. + /cid ccone 2 get def % cid は cone の 番号. + ccone 0 get /ccone set + { + ccone tag 0 eq { exit } { } ifelse + +% ccone の kk 番目の facet について flip する. + ccone kk cone.epsilon flipWeight /w set + (Trying new weight is ) messagen w message + w liftWeight /next_weight_w_wv set + (Trying new weight [w,wv] is ) messagen next_weight_w_wv message + + cone.ckmFlip { + [ + cone.gblist cid get (grobnerBasis) getNode 2 get % reduce gb + cone.vv + cone.gblist cid get (weight) getNode [2 0 2] get % weight + ccone kk getaVectorOnFacet liftWeight 1 get % weight on facet + next_weight_w_wv 1 get % new weight + ] /ttt set + ttt message + ttt ckmFlip /cone.cgb set + }{ + cone.input next_weight_w_wv 1 get cone.gb /cone.cgb set + } ifelse + + cone.cgb tag 0 eq not { + [w] next_weight_w_wv join /cone.cgb_weight set + next_weight_w_wv 1 get cone.cgb coneEq /cone.g_ineq set + cone.g_ineq cone.w_ineq join cone.Wt mul cone.Lpt mul + pruneZeroVector /cone.gw_ineq_projectedWtLpt set + + (cone.gw_ineq_projectedWtLpt is obtained.) message + + cone.gw_ineq_projectedWtLpt getConeInfo /cone.nextConeInfo set +% 次元を調べる. だめなら retry + cone.nextConeInfo 0 get 0 get to_int32 cone.d eq { + cone.nextConeInfo 1 get newCone /ncone set + ccone ncone getCommonFacet 0 get { + (Flip succeeded.) message + exit + } { } ifelse + } { } ifelse +% common face がなければ やはり epsilon を小さく. + cone.nextConeInfo 0 get 0 get to_int32 cone.d eq { + (ccone and ncone do not have a common facet.) message + } { + (ncone is not maximal dimensional. ) message + } ifelse + }{ } ifelse + + (Decreasing epsilon to ) messagen + cone.epsilon (1).. (2).. div mul /cone.epsilon set + cone.epsilon cone.epsilon.limit sub numerator (0).. lt { + (Too small cone.epsilon ) error + } { } ifelse + cone.epsilon message + } loop + /arg1 ncone def + ] pop + popVariables + arg1 +} def + +%%change +/cone_ir_input { + /arg1 set + [/msg ] pushVariables + [ + /msg arg1 def + (---------------) message + msg message + ( ) message + (Please also refer to the value of the variables cone.getConeInfo.rr0) message + ( cone.getConeInfo.rr1 cone.Lp cone.cinit) message + $ cone.cinit (FACETS) getNode :: $ message + (We are sorry that we cannot accept this input.) error + ] pop + popVariables } def