File: [local] / OpenXM / src / kan96xx / Doc / gfan.sm1 (download)
Revision 1.1, Sun Sep 5 10:19:29 2004 UTC (20 years ago) by takayama
Branch: MAIN
The program gfan.sm1 is for enumerating all maximal dimensional
Grobner cones of a given ideal. (a giant commit)
|
% $OpenXM: OpenXM/src/kan96xx/Doc/gfan.sm1,v 1.1 2004/09/05 10:19:29 takayama Exp $
% cp cone.sm1 $OpenXM_HOME/src/kan96xx/Doc/gfan.sm1
% $Id: cone.sm1,v 1.29 2004/09/05 10:16:00 taka Exp $
% iso-2022-jp
/cone.debug 1 def
/ox.k0.loaded boundp {
} {
[(parse) (ox.sm1) pushfile] extension
} ifelse
% Global: cone.type
% $B$I$N(B exponents $B$r<h$j=P$9$N$+;XDj$9$k(B.
% cf. exponents, gbext h $B$d(B H $B$b8+$k$+(B?
% 0 : x,y,Dx,Dy
% 1 : x,y,Dx,Dy,h,H
% 2 : x,y,Dx,Dy,h
/cone.type 2 def
% Global: cone.local
% cone.local: Local $B$+(B? 1 $B$J$i(B local
/cone.local 1 def
% Global: cone.h0
% cone.h0: 1 $B$J$i(B h $B$N(B weight 0 $B$G$N(B Grobner fan $B$r7W;;$9$k(B.
/cone.h0 1 def
% Global: cone.n (number of variables in GB)
% cone.m (freedom of the weight space. cf. cone.W)
% cone.d (pointed cones lies in this space. cf. cone.Lp)
% These are set during getting the cone.startingCone
%<
% Usage: wv g coneEq1
% in(f) $B$,(B monomial $B@lMQ(B. in_w(f) = LT(f) $B$H$J$k(B weight w $B$NK~$?$9(B
% $BITEy<0@)Ls$r5a$a$k(B.
%>
/coneEq1 {
/arg1 set
[/g /eqs /gsize /i /j /n /f /exps /m % Do not use "eq" as a variable
/expsTop
] pushVariables
[
/g arg1 def % Reduced Grobner basis
/eqs [ ] def % $BITEy<07O$N78?t(B
/gsize g length def
0 1 gsize 1 sub {
/i set
g i get /f set % f $B$O(B i $BHVL\$N(B reduced Grobner basis $B$N85(B
[(exponents) f cone.type] gbext /exps set % exps $B$O(B f $B$N(B exponent vector
exps length /m set
m 1 eq not {
/expsTop exps 0 get def % expsTop $B$O(B f $B$N@hF,$N(B exponent vector.
1 1 exps length 1 sub {
/j set
eqs [expsTop exps j get sub] join /eqs set
% exps[0]-exps[j] $B$r(B eqs $B$X3JG<$7$F$$$/$@$1(B.
% Cone $B$N(B closure $B$r$@$9$N$G(B >= $B$G(B OK.
} for
} { } ifelse
} for
/arg1 eqs def
] pop
popVariables
arg1
} def
%<
% Usage: ww g coneEq
% ww $B$O(B [v1 w1 v2 w2 ... ] $B7A<0(B. (v-w $B7A<0(B) w1, w2 $B$O(B univNumber $B$G$b$$$$(B.
% g $B$O(B reduced Grobner basis
% in(f) $B$,(B monomial $B$G$J$$>l9g$b07$&(B.
% in_w(f) = in_ww(f) $B$H$J$k(B weight w $B$NK~$?$9(B
% $BITEy<0@)Ls$r5a$a$k(B.
% ord_w, init (weightv) $B$rMQ$$$k(B.
%>
/coneEq {
/arg2 set
/arg1 set
[/g /eqs /gsize /i /j /n /f /exps /m
/expsTop /ww /ww2 /iterms
] pushVariables
[
/g arg2 def % Reduced Grobner basis
/ww arg1 def % weight vector. v-w $B7A<0(B
ww to_int32 /ww set % univNum $B$,$"$l$P(B int32 $B$KD>$7$F$*$/(B.
/ww2 ww weightv def % v-w $B7A<0$r(B $B?t;z$N%Y%/%H%k$K(B. (init $BMQ(B)
/eqs [ ] def % $BITEy<07O$N78?t(B
/gsize g length def
0 1 gsize 1 sub {
/i set
g i get /f set % f $B$O(B i $BHVL\$N(B reduced Grobner basis $B$N85(B
[(exponents) f cone.type] gbext /exps set % exps $B$O(B f $B$N(B exponent vector
exps length /m set
m 1 eq not {
/expsTop exps 0 get def % expsTop $B$O(B f $B$N@hF,$N(B exponent vector.
/iterms f ww2 init length def % f $B$N(B initial term $B$N9`$N?t(B.
% in_ww(f) > f_j $B$H$J$k9`$N=hM}(B.
iterms 1 exps length 1 sub {
/j set
eqs [expsTop exps j get sub] join /eqs set
% exps[0]-exps[j] $B$r(B eqs $B$X3JG<$7$F$$$/(B.
} for
% in_ww(f) = f_j $B$H$J$k9`$N=hM}(B.
[(exponents) f ww2 init cone.type] gbext /exps set % exps $B$O(B in(f)
1 1 iterms 1 sub {
/j set
eqs [exps j get expsTop sub] join /eqs set
eqs [expsTop exps j get sub] join /eqs set
% exps[j]-exps[0], exps[0]-exps[j] $B$r3JG<(B.
% $B7k2LE*$K(B (exps[j]-exps[0]).w = 0 $B$H$J$k(B.
} for
} { } ifelse
} for
/arg1 eqs def
] pop
popVariables
arg1
} def
%<
% Usage: wv g coneEq genPo
% polymake $B7A<0$N(B INEQUALITIES $B$r@8@.$9$k(B. coneEq -> genPo $B$HMxMQ(B
%>
/genPo {
/arg1 set
[/outConeEq /rr /nn /ii /mm /jj /ee] pushVariables
[
/outConeEq arg1 def
/rr [(INEQUALITIES) nl] cat def % $BJ8;zNs(B rr $B$KB-$7$F$$$/(B.
outConeEq length /nn set
0 1 nn 1 sub {
/ii set
outConeEq ii get /ee set
[ rr
(0 ) % $BHs$;$$$8MQ$N(B 0 $B$r2C$($k(B.
0 1 ee length 1 sub {
/jj set
ee jj get toString ( )
} for
nl
] cat /rr set
} for
/arg1 rr def
] pop
popVariables
arg1
} def
%<
% Usage: wv g coneEq genPo2
% doPolyamke $B7A<0$N(B INEQUALITIES $B$r@8@.$9$k(B. coneEq -> genPo2 $B$HMxMQ(B
% tfb $B7A<0J8;zNs(B.
%>
/genPo2 {
/arg1 set
[/outConeEq /rr /nn /ii /mm /jj /ee] pushVariables
[
/outConeEq arg1 def
/rr $polymake.data(polymake.INEQUALITIES([$ def
% $BJ8;zNs(B rr $B$KB-$7$F$$$/(B.
outConeEq length /nn set
0 1 nn 1 sub {
/ii set
outConeEq ii get /ee set
[ rr
([0,) % $BHs$;$$$8MQ$N(B 0 $B$r2C$($k(B.
0 1 ee length 1 sub {
/jj set
ee jj get toString
jj ee length 1 sub eq { } { (,) } ifelse
} for
(])
ii nn 1 sub eq { } { (,) } ifelse
] cat /rr set
} for
[rr $]))$ ] cat /rr set
/arg1 rr def
] pop
popVariables
arg1
} def
/test1 {
[(x,y) ring_of_differential_operators 0] define_ring
[ (x + y + Dx + Dy).
(x ^2 Dx^2 + y^2 Dy^2).
(x).
] /gg set
gg coneEq1 /ggc set
gg message
ggc pmat
ggc genPo message
} def
/test2 {
[(parse) (dhecart.sm1) pushfile] extension
dh.test.p1 /ff set
ff 0 get coneEq1 /ggc set
ggc message
ggc genPo /ss set
ss message
(Data is in ss) message
} def
/test3 {
% [(parse) (cohom.sm1) pushfile] extension
/ww [(Dx) 1 (Dy) 1] def
[(x,y) ring_of_differential_operators
[ww] weight_vector
0] define_ring
[ (x Dx + y Dy -1).
(y^2 Dy^2 + 2 + y Dy ).
] /gg set
gg {homogenize} map /gg set
[gg] groebner 0 get /gg set
ww message
ww gg coneEq /ggc set
gg message
ggc pmat
ggc genPo message
} def
%<
% Usage: test3b
% Grobner cone $B$r7hDj$7$F(B, polymake $BMQ$N%G!<%?$r@8@.$9$k%F%9%H(B.
% weight (0,0,1,1) $B$@$H(B max dim cone $B$G$J$$(B.
%>
/test3b {
% [(parse) (cohom.sm1) pushfile] extension
/ww [(Dx) 1 (Dy) 2] def
[(x,y) ring_of_differential_operators
[ww] weight_vector
0] define_ring
[ (x Dx + y Dy -1).
(y^2 Dy^2 + 2 + y Dy ).
] /gg set
gg {homogenize} map /gg set
[gg] groebner 0 get /gg set
ww message
ww gg coneEq /ggc set
gg message
ggc pmat
% ggc genPo /ggs set % INEQ $B$rJ8;zNs7A<0$G(B
% ggs message
% ggs output
% (mv sm1out.txt test3b.poly) system
% (Type in polymake-pear.sh test3b.poly FACETS) message
ggc genPo2 /ggs set % INEQ $B$rJ8;zNs7A<0(B for doPolymake
ggs message
} def
% commit (dr.sm1): lcm, denominator, ngcd, to_univNum, numerator, reduce
% 8/22, changelog-ja $B$^$@(B.
% to do : nnormalize_vec, sort_vec --> shell $B$G(B OK.
% 8/27, getNode
/test4 {
$polymake.data(polymake.INEQUALITIES([[0,1,0,0],[0,0,1,0]]))$ /ff set
[(FACETS) ff] doPolymake /rr set
rr 1 get /rr1 set
rr1 getLinearitySubspace pmat
} def
%<
% Usage: vv ineq isInLinearSpace
% vv $B$,(B ineq[i] > 0 $B$GDj5A$5$l$kH>6u4V$N$I$l$+$K$O$$$C$F$$$k$J$i(B 0
% vv $B$,(B $BA4$F$N(B i $B$K$D$$$F(B ineq[i] = 0 $B$K$O$$$C$F$$$?$i(B 1.
%>
/isInLinearSpace {
/arg2 set
/arg1 set
[/vv /ineq /ii /rr] pushVariables
[
/vv arg1 def
/ineq arg2 def
/rr 1 def
{
0 1 ineq length 1 sub {
/ii set
% vv . ineq[ii] != 0 $B$J$i(B vv $B$O(B linearity space $B$N85$G$J$$(B.
vv ineq ii get mul to_univNum isZero {
} { /rr 0 def exit} ifelse
} for
exit
} loop
/arg1 rr def
] pop
popVariables
arg1
} def
%<
% Usages: doPolymakeObj getLinearitySubspace
% INEQUALITIES $B$H(B VERTICES $B$+$i(B maximal linearity subspace
% $B$N@8@.%Y%/%H%k$r5a$a$k(B.
% $BNc(B: VERTICES [[0,1,0,0],[0,0,1,0],[0,0,0,-1],[0,0,0,1]]]
% $BNc(B: INEQUALITIES [[0,1,0,0],[0,0,1,0]]
% $BF~NO$O(B polymake $B$N(B tree (doPolymake $B$N(B 1 get)
%>
/getLinearitySubspace {
/arg1 set
[/pdata /vv /ineq /rr /ii] pushVariables
[
/pdata arg1 def
{
/rr [ ] def
% POINTED $B$J$i(B max lin subspace $B$O(B 0.
pdata (POINTED) getNode tag 0 eq { } { exit} ifelse
pdata (INEQUALITIES) getNode 2 get 0 get /ineq set
pdata (VERTICES) getNode 2 get 0 get /vv set
0 1 vv length 1 sub {
/ii set
% -vv[ii] $B$,(B ineq $B$rK~$?$9$+D4$Y$k(B.
vv ii get ineq isInLinearSpace {
rr [vv ii get] join /rr set
} { } ifelse
} for
exit
} loop
/arg1 rr def
] pop
popVariables
arg1
} def
%<
% Usages: mm asir_matrix_image
% $B@8@.85$h$j@~7A6u4V$N4pDl$rF@$k(B.
%>
/asir_matrix_image {
/arg1 set
[/mm /rr] pushVariables
[(CurrentRingp)] pushEnv
[
/mm arg1 def
mm to_univNum /mm set
oxasir.ccc [ ] eq {
(Starting ox_asir server.) message
ox_asirConnectMethod
} { } ifelse
{
oxasir.ccc [(matrix_image) mm] asir
/rr set
rr null_to_zero /rr set
exit
(asir_matrix_image: not implemented) error exit
} loop
rr numerator /rr set
/arg1 rr def
] pop
popEnv
popVariables
arg1
} def
[(asir_matrix_image)
[(Calling the function matrix_image of asir. It gets a reduced basis of a given matrix.)
(Example: [[1 2 3] [2 4 6]] asir_matrix_image)
]] putUsages
%<
% Usages: mm asir_matrix_kernel
% $BD>8r$9$k6u4V$N4pDl(B.
%>
/asir_matrix_kernel {
/arg1 set
[/mm /rr] pushVariables
[(CurrentRingp)] pushEnv
[
/mm arg1 def
mm to_univNum /mm set
oxasir.ccc [ ] eq {
(Starting ox_asir server.) message
ox_asirConnectMethod
} { } ifelse
{
oxasir.ccc [(matrix_kernel) mm] asir
/rr set
rr null_to_zero /rr set
exit
(asir_matrix_image: not implemented) error exit
} loop
rr 1 get numerator /rr set
/arg1 rr def
] pop
popEnv
popVariables
arg1
} def
[(asir_matrix_kernel)
[(Calling the function matrix_kernel of asir.)
(It gets a reduced basis of the kernel of a given matrix.)
(Example: [[1 2 3] [2 4 6]] asir_matrix_kernel)
]] putUsages
%<
% Usages: v null_to_zero
%>
/null_to_zero {
/arg1 set
[/pp /rr] pushVariables
[
/pp arg1 def
{
/rr pp def
pp isArray {
pp {null_to_zero} map /rr set
exit
}{ } ifelse
pp tag 0 eq {
/rr (0).. def
exit
}{ } ifelse
exit
} loop
/arg1 rr def
] pop
popVariables
arg1
} def
[(null_to_zero)
[(obj null_to_zero rob)
$It translates null to (0)..$
]] putUsages
% [2 0] lcm $B$O(B 0 $B$r$b$I$9$,$$$$$+(B? --> OK.
%<
% Usages: mm addZeroForPolymake
% $B0J2<$NFs$D$N4X?t$O(B, toQuotientSpace $B$K$bMxMQ(B.
% Polymake INEQUALITIES $BMQ$K(B 0 $B$r;O$a$KB-$9(B.
% $BF~NO$O(B $B%j%9%H$N%j%9%H(B
% [[1,2], [3,4],[5,6]] --> [[0,1,2],[0,3,4],[0,5,6]]
%>
/addZeroForPolymake {
/arg1 set
[/mm /rr] pushVariables
[
/mm arg1 def
mm to_univNum /mm set
mm { [(0)..] 2 1 roll join } map /mm set
/arg1 mm def
] pop
popVariables
arg1
} def
%<
% Usages: mm cone.appendZero
%>
/cone.appendZero {
/arg1 set
[/mm /rr] pushVariables
[
/mm arg1 def
mm to_univNum /mm set
mm { [(0)..] join } map /mm set
/arg1 mm def
] pop
popVariables
arg1
} def
%<
% Usages: mm removeFirstFromPolymake
% $B;O$a$N(B 0 $B$r<h$j=|$/(B.
% $BF~NO$O(B $B%j%9%H$N%j%9%H(B
% [[0,1,2],[0,3,4],[0,5,6]] ---> [[1,2], [3,4],[5,6]]
%>
/removeFirstFromPolymake {
/arg1 set
[/mm /rr] pushVariables
[
/mm arg1 def
mm to_univNum /mm set
mm {rest} map /mm set
/arg1 mm def
] pop
popVariables
arg1
} def
%<
% Usages: mm genUnit
% [1,0,0,...] $B$r2C$($k$?$a$K@8@.(B.
% [[0,1,2], [0,3,4],[0,5,6]]--> [1,0,0]
%>
/genUnit {
/arg1 set
[/mm /rr /i] pushVariables
[
/mm arg1 def
mm 0 get length newVector /rr set
rr null_to_zero /rr set
rr 0 (1).. put
/arg1 rr def
] pop
popVariables
arg1
} def
%<
% Usages: mm genUnitMatrix
% [[0,1,2], [0,3,4],[0,5,6]]--> [[1,0,0],[0,1,0],[0,0,1]]
%>
/genUnitMatrix {
/arg1 set
[/mm /rr /nn /i] pushVariables
[
/mm arg1 def
mm 0 get length /nn set
[
0 1 nn 1 sub {
/i set
nn newVector null_to_zero /mm set
mm i (1).. put
mm
} for
]
/arg1 set
] pop
popVariables
arg1
} def
%<
%%note: 2004, 8/29 (sun)
% toQuotientSpace : Linearity space $B$G3d$k(B.
% Usages: ineq mm toQuotientSpace
% $BF~NO$O(B coneEq $B$N=PNO(B ineq
% $B$*$h$S(B doPolymake --> getLinearitySubspace ==> L
% [L,[1,0,0,...]] asir_matrix_kernel removeFirstFromPolymake $B$GF@$i$l$?(B mm
% $B=PNO$+$i(B 0 $B%Y%/%H%k$O:o=|(B.
% $B=PNO$b(B coneEq $B7A<0(B. $BFC$K(B polymake $BMQ$K(B 0 $B$r2C$($k$N$,I,MW(B.
% ref: getUnit, removeFirstFromPolymake, addZeroForPolymake,
% asir_matrix_kernel, getLinearitySubspace
%>
/toQuotientSpace {
/arg2 set
/arg1 set
[/ineq /mm /rr] pushVariables
[
/ineq arg1 def
/mm arg2 def
ineq mm transpose mul /rr set
/arg1 rr def
] pop
popVariables
arg1
} def
/test5.data
$polymake.data(polymake.INEQUALITIES([[0,1,-1,1,-1,0],[0,0,-1,0,-1,2],[0,0,-1,0,-1,2],[0,0,-2,0,-2,4],[0,-1,0,-1,0,2],[0,-2,0,-2,0,4]]),polymake.VERTICES([[0,0,-1,0,0,0],[0,-1,-1,0,0,0],[0,1,0,-1,0,0],[0,-1,0,1,0,0],[0,0,1,0,-1,0],[0,0,-1,0,1,0],[0,-2,-2,0,0,-1],[0,2,2,0,0,1]]),polymake.FACETS([[0,1,-1,1,-1,0],[0,-1,0,-1,0,2]]),polymake.AFFINE_HULL(),polymake.FEASIBLE(),polymake.NOT__POINTED(),polymake.FAR_FACE([polymake._set([0,1,2,3,4,5,6,7])]),polymake.VERTICES_IN_INEQUALITIES([polymake._set([1,2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7])]),polymake.DIM([[5]]),polymake.AMBIENT_DIM([[5]]))$
def
%<
% Usages: test5
%% getConeInfo $B$rJQ99$9$l$P(B polymake $B$r8F$P$:$K%F%9%H$G$-$k(B.
%>
/test5 {
% test3b $B$h$j(B
/ww [(Dx) 1 (Dy) 2] def
% /ww [(x) 1 (y) -2 (Dx) 3 (Dy) 6] def
[(x,y) ring_of_differential_operators
[ww] weight_vector
0] define_ring
[ (x Dx + y Dy -1).
(y^2 Dy^2 + 2 + y Dy ).
] /gg set
gg {homogenize} map /gg set
[(AutoReduce) 1] system_variable
[gg] groebner 0 get /gg set
ww message
ww gg coneEq getConeInfo /rr set
(Type in rr 0 get :: ) message
} def
%[5, [[1,0,1,0,-2],[0,1,0,1,-2]], $NOT__POINTED$ ]
% $B$3$N>l9g$O(B 2 $B<!85$^$GMn$9$H(B pointed cone $B$K$J$k(B.
% coneEq mmc transpose $B$r$b$H$K(B FACETS $B$r7W;;$9$l$P$h$$(B.
%<
% Usage: ceq getConeInfo
% vw $B$O(B [v1 w1 v2 w2 ... ] $B7A<0(B. (v-w $B7A<0(B) w1, w2 $B$O(B univNumber $B$G$b$$$$(B.
% g $B$O(B reduced Grobner basis $B$H$7$F(B vw g coneEq $B$r7W;;(B. $B$3$l$r(B getConeInfo $B$X(B.
% Grobner cone $B$N(B $B<!85(B cdim (DIM), $BJd6u4V(B (linearity space ) $B$X$N9TNs(B mmc
% linearity space $B<+BN(B, pointed or not__pointed
% $B$D$^$j(B [cdim, L', L, PointedQ]
% $B$r7W;;$7$FLa$9(B. (polymake $B7A<0$NM>J,$JItJ,$J$7(B)
% polymake $BI,MW(B.
% ref: coneEq
% Global:
% cone.getConeInfo.rr0, cone.getConeInfo.rr1 $B$K(B polymake $B$h$j$NLa$jCM$,$O$$$k(B.
%>
/getConeInfo {
/arg1 set
[/ww /g /ceq /ceq2 /cdim /mmc /mmL /rr /ineq /ppt] pushVariables
[
/ceq arg1 def
ceq pruneZeroVector /ceq set
ceq genPo2 /ceq2 set
% ceq2 $B$O(B polymake.data(polymake.INEQUALITIES(...)) $B7A<0(B
% polymake $B$G(B ceq2 $B$N<!85$N7W;;(B.
/getConeInfo.ceq ceq def /getConeInfo.ceq2 ceq2 def
cone.debug { (Calling polymake DIM.) message } { } ifelse
[(DIM) ceq2] doPolymake 1 get /rr set
cone.debug {(Done.) message } { } ifelse
% test5 $B$K$O<!$N%3%a%s%H$H$j$5$k(B. $B>e$N9T$r%3%a%s%H%"%&%H(B.
% test5.data tfbToTree /rr set
/cone.getConeInfo.rr0 rr def
rr (DIM) getNode /cdim set
cdim 2 get 0 get 0 get 0 get to_univNum /cdim set
% polymake $B$N(B DIM $B$O0l$D>.$5$$$N$G(B 1 $BB-$9(B.
cdim (1).. add /cdim set
rr (FACETS) getNode tag 0 eq {
% FACETS $B$r;}$C$F$$$J$$$J$i:FEY7W;;$9$k(B.
% POINTED, NOT__POINTED $B$bF@$i$l$k(B
cone.debug { (Calling polymake FACETS.) message } { } ifelse
[(FACETS) ceq2] doPolymake 1 get /rr set
cone.debug { (Done.) message } { } ifelse
} { } ifelse
rr (VERTICES) getNode tag 0 eq {
(internal error: VERTICES is not found.) error
} { } ifelse
/cone.getConeInfo.rr1 rr def
rr (NOT__POINTED) getNode tag 0 eq {
% cone $B$,(B pointed $B$N;~$O(B mmc $B$OC10L9TNs(B. genUnitMatrix $B$r;H$&(B.
% VERTICES $B$h$j0l$D>.$5$$%5%$%:(B.
/mmc
[ rr (VERTICES) getNode 2 get 0 get 0 get rest]
genUnitMatrix
def
/mmL [ ] def
/ppt (POINTED) def
} {
% pointed $B$G$J$$>l9g(B,
% cone $B$N@~7AItJ,6u4V$r7W;;(B.
rr getLinearitySubspace /mmL set
[mmL genUnit] mmL join /mmc set % [1,0,0,...] $B$rB-$9(B.
mmc asir_matrix_kernel /mmc set % $BJd6u4V(B
mmc removeFirstFromPolymake /mmc set % $B$R$H$D>.$5$$%5%$%:$K(B.
[mmL genUnit] mmL join asir_matrix_image
removeFirstFromPolymake /mmL set
mmL asir_matrix_image /mmL set % Linearity space $B$r5a$a$k(B. rm 0vector
/ppt (NOT__POINTED) def
} ifelse
/arg1 [[cdim mmc mmL ppt] rr] def
] pop
popVariables
arg1
} def
/test.put {
/dog [(dog) [[(legs) 4] ] [1 2 3 ]] [(class) (tree)] dc def
/man [(man) [[(legs) 2] ] [1 2 3 ]] [(class) (tree)] dc def
/ma [(mammal) [ ] [man dog]] [(class) (tree)] dc def
/fan [ma 1 copy] def
ma (dog) getNode /dd set
dd 2 get /dd2 set
dd2 1 0 put
ma message
fan message
} def
/test6.data
$polymake.data(polymake.INEQUALITIES([[0,1,-1,1,-1,0],[0,0,-1,0,-1,2],[0,0,-1,0,-1,2],[0,0,-2,0,-2,4],[0,-1,0,-1,0,2],[0,-2,0,-2,0,4]]),polymake.VERTICES([[0,0,-1,0,0,0],[0,-1,-1,0,0,0],[0,1,0,-1,0,0],[0,-1,0,1,0,0],[0,0,1,0,-1,0],[0,0,-1,0,1,0],[0,-2,-2,0,0,-1],[0,2,2,0,0,1]]),polymake.FACETS([[0,1,-1,1,-1,0],[0,-1,0,-1,0,2]]),polymake.AFFINE_HULL(),polymake.FEASIBLE(),polymake.NOT__POINTED(),polymake.FAR_FACE([polymake._set([0,1,2,3,4,5,6,7])]),polymake.VERTICES_IN_INEQUALITIES([polymake._set([1,2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7]),polymake._set([0,2,3,4,5,6,7])]))$
def
% tfbToTree
/arrayToTree { [(class) (tree)] dc } def
%<
% polymake $B$h$jF@$i$l$?(B TreeObject $B$+$i(B TreeObject cone $B$r@8@.$9$k(B.
% Usages: test6.data tfbToTree newCone $B$GF0:n%F%9%H(B
%>
/test6 {
test6.data tfbToTree /rr set
rr newCone /rr2 set
} def
%<
% Usages: doPolymakeObj newCone
%>
/newCone {
/arg1 set
[/polydata /cone /facets /vertices /flipped /ineq
/facetsv /rr] pushVariables
[
/polydata arg1 def
polydata (FACETS) getNode tag 0 eq {
(newCone : no FACETS data.) error
} { } ifelse
% facets $B$OM-M}?t$N>l9g@55,2=$9$k(B. data/test11 $B$G(B $BM-M}?t$G$k(B.
polydata (FACETS) getNode 2 get 0 get to_univNum
{ nnormalize_vec} map /facets set
[[ ] ] facets join shell rest removeFirstFromPolymake /facets set
% vertices $B$O(B cone $B$N>e$K$"$k$N$G@0?tG\(B OK. $B@55,$+$9$k(B.
polydata (VERTICES) getNode 2 get 0 get to_univNum
{ nnormalize_vec} map /vertices set
[[ ] ] vertices join shell rest removeFirstFromPolymake /vertices set
% inequalities $B$OM-M}?t$N>l9g@55,2=$9$k(B.
polydata (INEQUALITIES) getNode 2 get 0 get to_univNum
{ nnormalize_vec } map /ineq set
[[ ] ] ineq join shell rest removeFirstFromPolymake /ineq set
[(cone) [ ]
[
[(facets) [ ] facets] arrayToTree
[(flipped) [ ] facets length newVector null_to_zero] arrayToTree
[(facetsv) [ ] facets vertices newCone_facetsv] arrayToTree
[(vertices) [ ] vertices] arrayToTree
[(inequalities) [ ] ineq] arrayToTree
]
] arrayToTree /cone set
/arg1 cone def
] pop
popVariables
arg1
} def
%<
% Usages: newCone_facetv
% facet vertices newCone_facetv
% facet $B$K$N$C$F$$$k(B vertices $B$r$9$Y$FNs5s(B.
%>
/newCone_facetv {
/arg2 set
/arg1 set
[/facet /vertices] pushVariables
[
/facet arg1 def /vertices arg2 def
[
0 1 vertices length 1 sub {
/ii set
facet vertices ii get mul isZero
{ vertices ii get } { } ifelse
} for
]
/arg1 set
] pop
popVariables
arg1
} def
%<
% Usages: newCone_facetsv
% facets vertices newCone_facetv
% facets $B$K$N$C$F$$$k(B vertices $B$r$9$Y$FNs5s(B. $B%j%9%H$r:n$k(B.
%>
/newCone_facetsv {
/arg2 set
/arg1 set
[/facets /vertices] pushVariables
[
/facets arg1 def /vertices arg2 def
facets { vertices newCone_facetv } map
/arg1 set
] pop
popVariables
arg1
} def
%<
% Usages: cone_random
%>
/cone_random.start (2).. def
/cone_random {
[(tdiv_qr)
cone_random.start (1103515245).. mul
(12345).. add
(2147483646)..
] mpzext 1 get /cone_random.start set
cone_random.start
} def
/cone_random.limit 40 def
/cone_random_vec {
/arg1 set
[/nn /rr] pushVariables
[
/nn arg1 def
[
0 1 nn 1 sub {
pop
[(tdiv_qr) cone_random cone_random.limit] mpzext 1 get
} for
] /arg1 set
] pop
popVariables
arg1
} def
%<
% Usages: getNewRandomWeight
%% max dim $B$N(B cone $B$r@8@.$9$k$?$a$K(B, random $B$J(B weight $B$r@8@.$9$k(B.
%% h, H $B$N=hM}$bI,MW(B.
%% $B@)Ls>r7o(B u+v >= 2t $B$r$_$?$9(B weight $B$,I,MW(B. $B$3$l$r$I$N$h$&$K:n$k$N$+(B?
%>
/getNewRandomWeight {
/arg1 set
[/vv /vvd /rr] pushVariables
[
/vv arg1 def
vv { (D) 2 1 roll 2 cat_n } map /vvd set
] pop
popVariables
arg1
} def
% test7 : univNum $B$N(B weight $B$,@5$7$/G'<1$5$l$k$+$N%F%9%H(B
% aux-cone.sm1
%<
% Usages: n d coneEqForSmallFan.2 (cone.type 2 $B@lMQ(B: x,y,Dx,Dy,h)
% n $BJQ?t$N?t(B, d zero $B$K$7$J$$JQ?t$N?t(B. d $B$O(B max dim cone $B$N<!85$H$J$k(B.
% $B$O$8$a$+$i(B d $B8D$NJQ?t(B.
% 4, 2 , s,t,x,y $B$J$i(B weight $B$O(B s,t,Ds,Dt $B$N$_(B.
% u_i + v_i >= 0 , u_i = v_i = 0.
% homog $BJQ?t$N>r7o(B u_i+v_i >= t, i.e, -t >= 0 $B$bF~$l$k(B.
% coneEq $B$N7k2L$H(B coneEqForSmallFan.2 $B$N7k2L$r(B join $B$7$F(B
% getConeInfo or newCone
% note-cone.sm1 2004.8.31 $B$r8+$h(B. w_ineq $B$"$?$j(B.
% cone.local $B$,@_Dj$5$l$F$$$k$H(B u_i <= 0 $B$b>r7o$KF~$k(B.
%>
/coneEqForSmallFan.2 {
/arg2 set
/arg1 set
[/n /d /nn /dd /ii /tt] pushVariables
[
/n arg1 def
/d arg2 def
n to_int32 /n set
d to_int32 /d set
/dd n d add def
/nn n n add def
% 0 ~ d-1, n ~ dd-1 $B$G$O(B u_i + v_i = 0
% d ~ n-1, dd ~ nn-1 $B$G$O(B u_i=v+i = 0.
% -t >= 0
[
% d ~ n-1, dd ~ nn-1 $B$G$O(B u_i=v+i = 0.
d 1 n 1 sub {
/ii set
% [ 0,0, ..., 0,1,0,... ; 0] $B$r@8@.(B
nn 1 add newVector null_to_zero /tt set
tt ii (1).. put
tt
% [ 0,0, ..., 0,-1,0,... ; 0] $B$r@8@.(B
nn 1 add newVector null_to_zero /tt set
tt ii (-1).. put
tt
} for
dd 1 nn 1 sub {
/ii set
nn 1 add newVector null_to_zero /tt set
tt ii (1).. put
tt
nn 1 add newVector null_to_zero /tt set
tt ii (-1).. put
tt
} for
% 0 ~ d-1, n ~ dd-1 $B$G$O(B u_i + v_i = 0
0 1 d 1 sub {
/ii set
nn 1 add newVector null_to_zero /tt set
tt ii (1).. put
tt ii n add (1).. put
tt
nn 1 add newVector null_to_zero /tt set
tt ii (-1).. put
tt ii n add (-1).. put
tt
} for
% -t >= 0
cone.h0 {
% t = 0
nn 1 add newVector null_to_zero /tt set
tt nn (1).. put
tt
nn 1 add newVector null_to_zero /tt set
tt nn (-1).. put
tt
}
{
% -t >= 0
nn 1 add newVector null_to_zero /tt set
tt nn (-1).. put
tt
} ifelse
% cone.local $B$,(B 1 $B$N;~(B
% 0 ~ d-1 $B$G$O(B -u_i >= 0
cone.local {
0 1 d 1 sub {
/ii set
nn 1 add newVector null_to_zero /tt set
tt ii (-1).. put
tt
} for
} { } ifelse
] /rr set
/arg1 rr to_univNum def
] pop
popVariables
arg1
} def
%<
% Usages: n d coneEqForSmallFan.1 (cone.type 1 $B@lMQ(B: x,y,Dx,Dy,h,H)
% cone.type 2 $B$G$O(B x,y,Dx,Dy,h
% coneEqForSmallFan.2 $B$N7k2L$rMQ$$$F@8@.(B.
% H $B$N>r7o$r2C$($k(B.
%>
/coneEqForSmallFan.1 {
/arg2 set
/arg1 set
[/n /d /i /j /rr /tt /tt2] pushVariables
[
/n arg1 def /d arg2 def
n d coneEqForSmallFan.2 /rr set
rr cone.appendZero /rr set
% H $BMQ$N(B 0 $B$r2C$($k(B.
% $B$H$j$"$($:(B t' = 0 $B$G$-$a$&$A(B.
cone.h0 { } { (cone.h0 = 0 has not yet been implemented.) error } ifelse
n 2 mul 2 add newVector null_to_zero /tt set
tt n 2 mul 2 add 1 sub (-1).. put
n 2 mul 2 add newVector null_to_zero /tt2 set
tt2 n 2 mul 2 add 1 sub (1).. put
rr [tt tt2] join /rr set
/arg1 rr to_univNum def
] pop
popVariables
arg1
} def
%<
% Usages: vv ineq toQuotientCone
% weight space $B$N(B $B%Q%i%a!<%?$D$1$N$?$a$K;H$&(B.
% cone.V $B$r5a$a$?$$(B. vv $B$O(B doPolymakeObj (VERTICES) getNode 2 get 0 get $B$GF@$k(B.
% vertices $B$N(B non-negative combination $B$,(B cone.
% vertice cone.w_ineq isInLinearSubspace $B$J$i<h$j=|$/(B.
% $B$D$^$j(B vertice*cone.w_ineq = 0 $B$J$i<h$j=|$/(B.
%
% $B$3$l$G@5$7$$(B? $B>ZL@$O(B? $B$^$@ESCf(B. cone.W $B$r5a$a$k$N$K;H$&(B. (BUG)
% cone.w_cone 1 get (VERTICES) getNode :: $B$HHf3S$;$h(B.
% $B$3$N4X?t$r8F$s$G(B cone.W $B$r:n$k$N$OITMW$+$b(B.
%
% Example: cf. parametrizeSmallFan
% 4 2 coneEqForSmallFan.2 /cone.w_ineq set cone.w_ineq getConeInfo /rr set
% rr 1 get (VERTICES) getNode 2 get 0 get removeFirstFromPolymake /vv set
% vv cone.w_ineq toQuotientCone pmat
%>
/toQuotientCone {
/arg2 set /arg1 set
[/vv /ineq /rr] pushVariables
[
/vv arg1 def /ineq arg2 def
vv {
dup
ineq isInLinearSpace 1 eq { pop }
{ } ifelse
} map /arg1 set
] pop
popVariables
arg1
} def
%<
% Usages: n d parametrizeSmallFan
% n : x $BJQ?t$N?t(B.
% d : 0 $B$K$7$J$$(B weight $B$N?t(B.
% $B<!$NBg0hJQ?t$b@_Dj$5$l$k(B.
% cone.W : weight $B$r%Q%i%a!<%?$E$1$9$k%Y%/%H%k$NAH(B.
% cone.Wpos : i $B$,(B 0 ~ Wpos-1 $B$NHO0O$N$H$-(B V[i] $B$X$O(B N $B$N85$r3]$1;;$7$F$h$$(B,
% i $B$,(B Wpos ~ $B$NHO0O$N$H$-(B V[i] $B$X$O(B Z $B$N85$r3]$1;;$7$F$h$$(B.
% cone.w_ineq : weight space $B$NITEy<0@)Ls(B. $B0J8e$N7W;;$G>o$KIU2C$9$k(B.
% cone.w_cone : w_ineq $B$r(B polymake $B$G(B getConeInfo $B$7$?7k2L(B.
% Example: /cone.local 1 def ; 4 2 parametrizeSmallFan pmat
% Example: /cone.local 0 def ; 4 2 parametrizeSmallFan pmat
%>
/parametrizeSmallFan {
/arg2 set /arg1 set
[/n /d /vv /coneray] pushVariables
[
/n arg1 def /d arg2 def
{
cone.type 1 eq {
n d coneEqForSmallFan.1 /cone.w_ineq set
exit
} { } ifelse
cone.type 2 eq {
n d coneEqForSmallFan.2 /cone.w_ineq set
exit
} { } ifelse
(This cone.type has not yet been implemented.) error
} loop
cone.w_ineq getConeInfo /cone.w_cone set
cone.w_cone 1 get (VERTICES) getNode 2 get 0 get
removeFirstFromPolymake /vv set
vv cone.w_ineq toQuotientCone /coneray set
coneray length /cone.Wpos set
coneray cone.w_cone 0 get 2 get join /cone.W set
/arg1 cone.W def
] pop
popVariables
arg1
} def
%<
% Usages: n d coneEqForTotalFan.2 (cone.type 2 $B@lMQ(B: x,y,Dx,Dy,h)
% n $BJQ?t$N?t(B,
% d 0 $B$K$7$J$$JQ?t(B.
% u_i + v_i >= 0 ,
% homog $BJQ?t$N>r7o(B u_i+v_i >= 0, t = 0 $B$bF~$l$k(B.
% coneEq $B$N7k2L$H(B coneEqForSmallFan.2 $B$N7k2L$r(B join $B$7$F(B
% getConeInfo or newCone
% cone.local $B$,@_Dj$5$l$F$$$k$H(B u_i <= 0 $B$b>r7o$KF~$k(B.
%>
/coneEqForTotalFan.2 {
/arg2 set
/arg1 set
[/n /nn /dd /ii /tt] pushVariables
[
/n arg1 def
/d arg2 def
n to_int32 /n set
d to_int32 /d set
/nn n n add def
/dd n d add def
% 0 ~ d-1, n ~ dd-1 $B$G$O(B u_i + v_i >= 0
% d ~ n-1, dd ~ nn-1 $B$G$O(B u_i=v+i = 0.
% t = 0
[
% d ~ n-1, dd ~ nn-1 $B$G$O(B u_i=v+i = 0.
d 1 n 1 sub {
/ii set
% [ 0,0, ..., 0,1,0,... ; 0] $B$r@8@.(B
nn 1 add newVector null_to_zero /tt set
tt ii (1).. put
tt
% [ 0,0, ..., 0,-1,0,... ; 0] $B$r@8@.(B
nn 1 add newVector null_to_zero /tt set
tt ii (-1).. put
tt
} for
dd 1 nn 1 sub {
/ii set
nn 1 add newVector null_to_zero /tt set
tt ii (1).. put
tt
nn 1 add newVector null_to_zero /tt set
tt ii (-1).. put
tt
} for
% 0 ~ d-1, n ~ dd-1 $B$G$O(B u_i + v_i >= 0
0 1 d 1 sub {
/ii set
nn 1 add newVector null_to_zero /tt set
tt ii (1).. put
tt ii n add (1).. put
tt
} for
% t = 0
cone.h0 {
% t = 0
nn 1 add newVector null_to_zero /tt set
tt nn (1).. put
tt
nn 1 add newVector null_to_zero /tt set
tt nn (-1).. put
tt
}
{
(coneForTotalFan.2. Not implemented.) error
} ifelse
% cone.local $B$,(B 1 $B$N;~(B
% 0 ~ d-1 $B$G$O(B -u_i >= 0
cone.local {
0 1 d 1 sub {
/ii set
nn 1 add newVector null_to_zero /tt set
tt ii (-1).. put
tt
} for
} { } ifelse
] /rr set
/arg1 rr to_univNum def
] pop
popVariables
arg1
} def
%<
% Usages: n d parametrizeTotalFan
% n : x $BJQ?t$N?t(B.
% d : 0 $B$K$7$J$$?t(B.
% $B<!$NBg0hJQ?t$b@_Dj$5$l$k(B.
% cone.W : weight $B$r%Q%i%a!<%?$E$1$9$k%Y%/%H%k$NAH(B.
% cone.Wpos : i $B$,(B 0 ~ Wpos-1 $B$NHO0O$N$H$-(B V[i] $B$X$O(B N $B$N85$r3]$1;;$7$F$h$$(B,
% i $B$,(B Wpos ~ $B$NHO0O$N$H$-(B V[i] $B$X$O(B Z $B$N85$r3]$1;;$7$F$h$$(B.
% cone.w_ineq : weight space $B$NITEy<0@)Ls(B. $B0J8e$N7W;;$G>o$KIU2C$9$k(B.
% cone.w_ineq $B$r(B getConeInfo $B$7$?7k2L$O(B cone.w_cone
% Example: /cone.local 1 def ; 3 parametrizeSmallFan pmat
% Example: /cone.local 0 def ; 3 parametrizeSmallFan pmat
% local $B$,(B 1 $B$@$H(B u_i <= 0 $B$K$J$k(B.
%>
/parametrizeTotalFan {
/arg2 set
/arg1 set
[/n /d /vv /coneray] pushVariables
[
/n arg1 def /d arg2 def
{
cone.type 2 eq { n d coneEqForTotalFan.2 /cone.w_ineq set exit}
{ } ifelse
(This cone.type has not yet been implemented.) error
} loop
cone.w_ineq getConeInfo /cone.w_cone set
cone.w_cone 1 get (VERTICES) getNode 2 get 0 get
removeFirstFromPolymake /vv set
vv cone.w_ineq toQuotientCone /coneray set
coneray length /cone.Wpos set
coneray cone.w_cone 0 get 2 get join /cone.W set
/arg1 cone.W def
] pop
popVariables
arg1
} def
%<
% Usages: vlist wlist cone_wtowv
% [x y Dx Dy h] [-1 0 1 0 0] ==> [(x) -1 (Dx) 1] $B$r:n$k(B.
%>
/cone_wtowv {
/arg2 set /arg1 set
[/vlist /wlist /ii] pushVariables
[
/vlist arg1 def
/wlist arg2 def
wlist length vlist length eq {
} { (cone_wtowv: length of the argument must be the same.) error} ifelse
wlist to_int32 /wlist set
[
0 1 wlist length 1 sub {
/ii set
wlist ii get 0 eq { }
{ vlist ii get wlist ii get } ifelse
} for
] /arg1 set
] pop
popVariables
arg1
} def
%<
% Usages: pruneZeroVector
% genPo, getConeInfo $BEy$NA0$K;H$&(B. 0 $B%Y%/%H%k$O0UL#$N$J$$@)Ls$J$N$G=|$/(B.
%>
/pruneZeroVector {
/arg1 set
[/mm /ii /jj /tt] pushVariables
[
/mm arg1 def
mm to_univNum /mm set
[
0 1 mm length 1 sub {
/ii set
mm ii get /tt set
{
0 1 tt length 1 sub {
/jj set
tt jj get (0).. eq { }
{ tt exit } ifelse
} for
exit
} loop
} for
] /arg1 set
] pop
arg1
} def
%<
% Usages: a projectIneq v , dim(a) = n, dim(v) = d
% a*cone.Wt*cone.Lpt
%>
/projectIneq {
cone.Wt mul cone.Lpt mul
} def
%<
% Usages: v liftWeight [w vw], dim(v) = d, dim(w) = n, vw : vw $B7A<0$N(B weight
% v*cone.Lp*cone.W cone.vlist w cone_wtowv
%>
/liftWeight {
/arg1 set
[/v /w /vw] pushVariables
[
/v arg1 def
v cone.Lp mul cone.W mul /w set
[w cone.vlist w cone_wtowv] /arg1 set
] pop
popVariables
arg1
} def
%<
% Usage: m isZero
% dr.sm1 $B$X0\$9(B.
%>
/isZero {
/arg1 set
[/mm /ans /ii] pushVariables
[
/mm arg1 def
/ans 1 def
mm isArray {
0 1 mm length 1 sub {
/ii set
mm ii get isZero /ans set
ans 0 eq { exit } { } ifelse
} for
} {
{
mm tag 1 eq {/ans mm 0 eq def exit} { } ifelse
mm isPolynomial { /ans mm (0). eq def exit } { } ifelse
mm isUniversalNumber { /ans mm (0).. eq def exit } { } ifelse
/ans 0 def exit
} loop
} ifelse
/arg1 ans def
] pop
popVariables
arg1
} def
[(isZero)
[(m isZero bool)]] putUsages
%<
% Usage: m isNonNegative
% dr.sm1 $B$X0\$9(B.
%>
/isNonNegative {
/arg1 set
[/mm /ans /ii] pushVariables
[
/mm arg1 def
/ans 1 def
mm isArray {
0 1 mm length 1 sub {
/ii set
mm ii get isNonNegative /ans set
ans 0 eq { exit } { } ifelse
} for
} {
{
mm tag 1 eq {/ans mm 0 gt mm 0 eq or def exit} { } ifelse
mm isUniversalNumber { /ans mm (0).. gt mm (0).. eq or def exit }
{ } ifelse
mm isRational { mm (numerator) dc mm (denominator) dc mul /mm set
/ans mm (0).. gt mm (0).. eq or def exit } { } ifelse
/ans 0 def exit
} loop
} ifelse
/arg1 ans def
] pop
popVariables
arg1
} def
[(isNonNegative)
[(m isNonNegative bool)
(In case of matrix, m[i,j] >= 0 must hold for all i,j.)
]] putUsages
% Global variable: cone.weightBorder
% /cone.weightBorder null def $BITMW$G$"$m$&(B. getStartingCone $B$G@_Dj$5$l$k(B.
%<
% Usages: cone i isOnWeigthBorder
% cone $B$N(B i $BHVL\$N(B facet $B$,(B weight $B6u4V$N6-3&$K$"$k$+(B?
% $BBg0hJQ?t(B cone.weightBorder $B$,@_Dj$5$l$F$k$3$H(B.
% $B$3$NJQ?t$O(B cone $B$N(B facet $B%Y%/%H%k$N%j%9%H(B.
% $B$3$NJQ?t$O(B setWeightBorder $B$G@_Dj(B
% cone.weightBorder[0] or cone.weightBorder[1] or ...
% /ccone cone.startingCone def ccone 0 isOnWeightBorder
% ccone 1 isOnWeightBorder
%>
/isOnWeightBorder {
/arg2 set /arg1 set
[/cone /facet_i /i /j /vv /co /ans] pushVariables
[
/cone arg1 def /facet_i arg2 def
facet_i to_int32 /facet_i set
/ans 0 def
cone (facetsv) getNode 2 get facet_i get /vv set % Facet $B$r(B vertex $BI=8=(B.
{
0 1 cone.weightBorder length 1 sub {
/i set
cone.weightBorder i get /co set % co $B$K@)Ls>r7o(B
vv cone.Lp mul % vv $B$r(B weight space $B$X(B lift.
co mul isZero
{ /ans 1 def exit } { } ifelse
} for
exit
} loop
/arg1 ans def
] pop
popVariables
arg1
} def
%<
% Usages: cone i markFlipped
% cone $B$N(B i $BHVL\$N(B facet $B$K(B flipped $B$N0u$r$D$1$k(B. cone $B<+BN$,JQ99$5$l$k(B.
% cone $B$O(B class-tree. Constructor $B$O(B newCone
%>
/markFlipped {
/arg2 set /arg1 set
[/cone /facet_i /vv] pushVariables
[
/cone arg1 def /facet_i arg2 def
facet_i to_int32 /facet_i set
cone (flipped) getNode 2 get /vv set
vv facet_i (1).. put
] pop
popVariables
} def
%<
% Usages: cone getNextFacet i
% flipped $B$N(B mark $B$N$J$$(B facet $B$N(B index facet_i $B$rLa$9(B.
% $B$=$l$,$J$$$H$-$O(B null
%>
/getNextFacet {
/arg1 set
[/cone /facet_i /vv /ii] pushVariables
[
/cone arg1 def
/facet_i null def
cone (flipped) getNode 2 get /vv set
0 1 vv length 1 sub {
/ii set
vv ii get to_int32 0 eq { /facet_i ii def exit }
{ } ifelse
} for
/arg1 facet_i def
] pop
popVariables
arg1
} def
%<
% Usages: cone i epsilon flipWeight
% cone $B$N(B i $BHVL\$N(B facet $B$K$+$s$7$F(B flip $B$9$k(B.
% $B?7$7$$(B weight $B$r5a$a$k(B. cf. liftWeight
%>
/flipWeight {
/arg3 set /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
/ep arg3 def
ep to_univNum (1).. div /ep set
% note: 2004.9.2
cone (facetsv) getNode 2 get facet_i get /v set
cone (facets) getNode 2 get facet_i get /f set
/vp v 0 get def
1 1 v length 1 sub {
/ii set
vp v ii get add /vp set
} for
vp ep f mul sub /vp set
vp nnormalize_vec /vp set
/arg1 vp def
] pop
popVariables
arg1
} def
%<
% Usages: cone1 cone2 isSameCone bool
% cone1 cone2 $B$,Ey$7$$$+(B? facet $B$GHf$Y$k(B.
% cone1, cone2 $B$O(B pointed cone $B$G$J$$$H$$$1$J$$(B.
%>
/isSameCone {
/arg2 set /arg1 set
[/cone1 /cone2 /facets1 /facets2 /ans] pushVariables
[
/cone1 arg1 def
/cone2 arg2 def
/facets1 cone1 (facets) getNode 2 get def
/facets2 cone2 (facets) getNode 2 get def
facets1 length facets2 length eq {
facets1 facets2 sub isZero /ans set
} {
/ans 0 def
} ifelse
/arg1 ans def
] pop
popVariables
arg1
} def
%<
% Usages: cone1 cone2 getCommonFacet list
% cone1 $B$NCf$G(B cone2 $B$K4^$^$l$k(B facet $B$N%j%9%H(B
% cone2 $B$NCf$G(B cone1 $B$K4^$^$l$k(B facet $B$N%j%9%H$r$b$I$9(B.
% [1 [i] [j]] $B$"$k$H$-(B. [0 [ ] [ ]] $B$J$$$H$-(B.
% cone1 $B$N(B facetsv[i] $B$,(B cone2 $B$K4^$^$l$k$+D4$Y$k(B.
% cone2 $B$N(B facetsv[i] $B$,(B cone1 $B$K4^$^$l$k$+D4$Y$k(B.
% cone1, cone2 $B$O(B pointed cone $B$G$J$$$H$$$1$J$$(B.
%>
/getCommonFacet {
/arg2 set /arg1 set
[/cone1 /cone2 /facets /ineq /ans1 /ans2 /i /tt] pushVariables
[
/cone1 arg1 def
/cone2 arg2 def
/facets cone1 (facetsv) getNode 2 get def
/ineq cone2 (inequalities) getNode 2 get def
/ans1 [
0 1 facets length 1 sub {
/i set
facets i get /tt set % facetsv[i] $B$r(B tt $B$X(B.
ineq tt transpose mul isNonNegative {
i
} { } ifelse
} for
] def
/facets cone2 (facetsv) getNode 2 get def
/ineq cone1 (inequalities) getNode 2 get def
/ans2 [
0 1 facets length 1 sub {
/i set
facets i get /tt set % facetsv[i] $B$r(B tt $B$X(B.
ineq tt transpose mul isNonNegative {
i
} { } ifelse
} for
] def
ans1 length 1 gt ans2 length 1 gt or {
(getCommonFacet found more than 1 common facets.) error
} { } ifelse
% $B6&DL(B facet $B$,$"$l$P(B 1, $B$J$1$l$P(B 0.
ans1 length 1 eq ans2 length 1 eq and {
/tt 1 def
} {
/tt 0 def
} ifelse
/arg1 [tt ans1 ans2] def
] pop
popVariables
arg1
} def
%
% -------------------------------------------------
% test8 $B$O(B aux-cone.sm1 $B$X0\F0(B.
% $B0J2<$$$h$$$h0lHL$N%W%m%0%i%`$N:n@.3+;O(B.
% -------------------------------------------------
%
%<
% Usages: setWeightBorder
% cone.weightBorder (weight cone $B$N(B facet $B%Y%/%H%k$N=89g(B) $B$r@_Dj$9$k(B.
% $B$"$HI{;:J*$H$7$F(B cone.w_cone_projectedWt (doPolymakeObj)
% cone.w_ineq_projectedWt
% cone.m $B<!85$N%Y%/%H%k(B.
% cone.W, cone.Wt, cone.w_ineq $B$,$9$G$K7W;;$:$_$G$J$$$H$$$1$J$$(B.
%>
/setWeightBorder {
[
(Entering setWeightBorder ) message
cone.w_ineq cone.Wt mul pruneZeroVector /cone.w_ineq_projectedWt set
{
cone.w_ineq_projectedWt length 0 eq {
% weight $B$N6u4V$K(B border $B$,$J$$>l9g(B.
/cone.weightBorder [ ] def
exit
} { } ifelse
% weight $B$N6u4V$K(B border $B$,$"$k>l9g(B.
cone.w_ineq_projectedWt getConeInfo /cone.w_cone_projectedWt set
cone.w_cone_projectedWt 0 get 0 get to_int32 cone.m to_int32 eq {
} {
(setWeightBorder : internal error.) message
} ifelse
cone.w_cone_projectedWt 1 get (FACETS) getNode 2 get 0 get
removeFirstFromPolymake /cone.weightBorder set
exit
} loop
(cone.weightBorder=) message
cone.weightBorder pmat
] pop
} def
%
% -------------------------------------------------
% $B%W%m%0%i%`$NN.$l(B.
% Global: cone.fan cone $B$rG[Ns$H$7$F3JG<$9$k(B.
%
% ncone (next cone) $B$,?75,$KF@$i$l$?(B cone $B$G$"$k$H$9$k(B.
% $B$3$N$H$-<!$NA`:n$r$9$k(B.
% 0. ncone $B$,(B cone.fan $B$K$9$G$K$J$$$+D4$Y$k(B. $B$"$l$P(B, internal error.
% 1. ncone markBorder ; ncone $B$NCf$N(B border $B>e$N(B facet $B$r(B mark
% 2. cone.fan $B$NCf$N(B cone $B$H6&DL(B facet $B$,$J$$$+D4$Y(B (getCommonFacet),
% $B$"$l$P$=$l$i$r(B mark $B$9$k(B.
% global: cone.incidence $B$K(B $B6&DL(Bfacet $B$r;}$DAH$_$N>pJs$r2C$($k(B.
% 3. ncone $B$r(B cone.fan $B$N:G8e$K2C$($k(B.
% $B0J>e$NA`:n$r$^$H$a$?$b$N$,(B ncone updateFan
%
% getNextFlip $B$O(B cone.fan $B$NCf$+$i(B flip $B$7$F$J$$(B cone $B$H(B facet $B$NAH$rLa$9(B.
% $B$J$1$l$P(B null $B$rLa$9(B. null $B$,La$l$P%W%m%0%i%`=*N;(B.
%
% getStargingCone $B$O7W;;$r=PH/$9$Y$-?75,$N(B cone $B$r7W;;$9$k(B. $BBg0hJQ?t(B cone.Lt, cone.W
% $B$J$I$b$3$NCf$G@_Dj$9$k(B.
% $BJQ?t%j%9%H(B, weight space $B$r@8@.$9$k4X?t(B, $BF~NOB?9`<0(B, weight $B$N8uJd(B $BEy$OBg0hJQ?t(B
% $B$H$7$FF~NO$7$F$*$/(B.
%
% reduced gb $B$O(B $B4X?t(B input weight cone.gb reduced_G $B$G7W;;$9$k(B.
%
%
% [ccone i] getNextCone ncone : flip $B$K$h$j<!$N(B cone $B$rF@$k(B.
%
% 1. clearGlobals ; $BF~NOBg0hJQ?t$N@_Dj(B.
% 2. getStartingCone /ncone set
% 3. { ncone updateFan
% 4. getNextFlip /cone.nextflip set
% 6. cone.nextflip isNull { exit } { } ifelse
% 7. cone.nextflip getNextCone /ncone set
% 8. } loop
%
%
% -------------------------------------------------
%
%<
% Usages: input weight cone.gb_Dh reduced_G
% gb in h[1,1](D)
%>
/cone.gb_Dh {
/arg2 set /arg1 set
[/ff /ww /gg] 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] groebner 0 get /gg set
/cone.gb_Dh.g gg def
/arg1 gg def
] pop
popVariables
arg1
} def
%<
% Usages: cone.boundp
%
/cone.boundp {
dup boundp 2 1 roll tag 0 eq not and
} def
%<
% Usages: clearGlobals
% cf. cone.boundp
% polymake $B$r:FEY8F$V$?$a$K(B global $BJQ?t$r%/%j%"$9$k(B.
% $B$^$@ESCf(B.
%>
/clearGlobals {
/cone.W null def
/cone.Wt null def
/cone.cinit null def
/cone.weightBorder null def
} def
%<
% Usages: getStartingCone ncone
% getStargingCone $B$O7W;;$r=PH/$9$Y$-?75,$N(B cone $B$r7W;;$9$k(B.
% $B@_Dj$9$Y$-Bg0hJQ?t$O0J2<$r8+$h(B.
%>
/getStartingCone.test {
%------------------Globals----------------------------------------
% --------------- $BF~NO%G!<%?MQBg0hJQ?t$N@_Dj(B --------------------------
%
% cone.input : $BF~NOB?9`<07O(B
/cone.input
[(t1-x-y) (h*t2-x^2-y^2) (2*x*Dt2+h*Dt1+h*Dx) (2*y*Dt2+h*Dt1+h*Dy)]
def
% cone.vlist : $BA4JQ?t$N%j%9%H(B
/cone.vlist [(t1) (t2) (x) (y) (Dt1) (Dt2) (Dx) (Dy) (h)] def
% cone.vv : define_ring $B7A<0$NJQ?t%j%9%H(B.
% t1,t2, x,y : t-space $B$N(B Grobner fan (local) $B$r5a$a$k(B.
/cone.vv (t1,t2,x,y) def
% cone.parametrizeWeightSpace : weight $B6u4V$r(B parametrize $B$9$k4X?t(B.
% $BBg0hJQ?t(B cone.W , cone.Wpos $B$b$-$^$k(B.
/cone.parametrizeWeightSpace {
4 2 parametrizeSmallFan
} def
% cone.w_start : weight$B6u4V$K$*$1$k(B weight $B$N=i4|CM(B.
% $B$3$NCM$G(B max dim cone $B$,F@$i$l$J$$$H(B random weight $B$K$h$k(B $B%5!<%A$,;O$^$k(B.
/cone.w_start
[ 1 4 ]
def
% cone.gb : gb $B$r7W;;$9$k4X?t(B.
/cone.gb {
cone.gb_Dh
} def
%
% ----------------- $B$*$o$j(B ---------------------------
%
} def % end of getStartingCone.test
/getStartingCone {
[/wv_start /w_start /reduced_G] pushVariables
[
% cone.n $B$O<+F0E*$K$-$a$i$l$k(B.
% cone.n $B$O(B GB $B$r7W;;$9$k6u4V$N<!85(B.
/cone.n cone.vlist length def
%[1] cone.W, cone.Wpos $B$r5a$a$k(B. cone.m $B$O(B cone.W $B$h$j<+F0E*$K$-$^$k(B.
% cone.m $B$O(B weight $B6u4V$N<+M3EY(B. cone.W $B$G<M1F$5$l$k@h$N<!85(B.
/cone.W cone.boundp {
(Skip cone.parametrizeWeightSpace. cf. clearGlobals) message
} {
cone.parametrizeWeightSpace
} ifelse
(parametrizing weight space: cone.W = ) messagen cone.W message
/cone.Wt cone.W transpose def
/cone.m cone.W length def
% WeightBorder $B$N>r7oH=Dj(B facet $B$r@_Dj(B.
/cone.weightBorder cone.boundp {
(Skip setWeightBorder cf. clearGlobals) message
} {
setWeightBorder
} ifelse
%[2] weight vector wv_start $B$r@8@.$9$k(B.
% wv_start $B$r@_Dj(B.
cone.w_start tag 0 eq {
% cone.w_start $B$,(B null $B$J$i(B random $B$K(B weight $B$r@_Dj(B.
/cone.w_start cone.m cone_random_vec def
} {
cone.w_start length cone.m to_int32 eq {
} {
(Error: cone.w_start has wrong length.) error
/cone.w_start cone.m cone_random_vec def
} ifelse
} ifelse
/w_start cone.w_start cone.W mul def
{
cone.vlist w_start cone_wtowv /wv_start set
(Trying a starting weight vector : ) messagen
wv_start pmat
%[3] reduced GB $B$N7W;;(B.
cone.input wv_start cone.gb /reduced_G set
(Reduced GB : ) message
reduced_G pmat
%[4] $B<M1F$7$F$+$i(B polytope $B$N%G!<%?$r7W;;(B.
wv_start reduced_G coneEq /cone.g_ineq set
cone.g_ineq cone.w_ineq join /cone.gw_ineq set
cone.gw_ineq cone.Wt mul /cone.gw_ineq_projectedWt set % $B<M1F(B
/cone.cinit cone.boundp {
(Skipping cone.gw_ineq_projectedWt getConeInfo. cf. clearGlobals) message
} {
cone.gw_ineq_projectedWt getConeInfo /cone.cinit set
} ifelse
(cone.cinit is --- the first number is the dim of cone.) messagen
cone.cinit 0 get pmat
% Maximal dimensional cone $B$+$I$&$+$N8!::(B. $B8!::$K%Q%9$9$l$P(B loop $B$r(B exit
% $B%Q%9$7$J$$>l9g(B w_start $B$r(B cone_random_vec $B$rMQ$$$FJQ99$9$k(B.
cone.cinit 0 get 0 get to_int32 cone.m eq { exit }
{
(Failed to get the max dim cone. Updating the weight ...) messagen
/w_start cone.m cone_random_vec cone.W mul def
% cone.cinit $B$r:FEY7W;;$9$k$?$a$K(B clear $B$9$k(B.
/cone.cinit null def
} ifelse
} loop
(cone.m = ) messagen cone.m message
(Suceeded to get the maximal dimensional startingCone.) message
% Linearity subspace $B$N(B orth complement $B$X$N<M1F9TNs(B.
% $BBg0hJQ?t(B cone.Lp, cone.Lpt $B$r@_Dj(B
cone.cinit 0 get 1 get /cone.Lp set
cone.Lp transpose /cone.Lpt set
% Linearity subspace $B$N9TNs$r@_Dj(B.
% $BBg0hJQ?t(B cone.L $B$r@_Dj(B
cone.cinit 0 get 2 get /cone.L set
% cone.d $B$O(B cone.W $B$*$h$S(B Linearity space $B$G3d$C$?8e(B, cone $B$r9M$($k$H$-$N<!85(B.
% $BBg0hJQ?t(B cone.d $B$N@_Dj(B.
/cone.d cone.Lp length def
cone.m cone.d eq {
(There is no linearity space) message
} {
(Dim of the linearity space is ) messagen cone.m cone.d sub message
(cone.Lp = ) messagen cone.Lp pmat
} ifelse
%[5] cone.g_ineq * cone.Wt * cone.Lpt
% cone.w_ineq * cone.Wt * cone.Lpt
% $B$G@)Ls$r(B d $B<!85%Y%/%H%k$KJQ49(B.
% W (R^m) $B6u4V$NITEy<0@)Ls$r(B L' (R^d) $B6u4V$X<M1F(B
% cone.gw_ineq_projectedWtLpt
% = cone.g_ineq*cone.Wt*cone.Lpt \/ cone.w_ineq*coneWt*cone.Lpt
/cone.gw_ineq_projectedWtLpt
cone.gw_ineq_projectedWt cone.Lpt mul
def
cone.m cone.d eq {
/cone.cinit.d cone.cinit def
} {
% cone.m > cone.d $B$J$i$P(B, $B:FEY(B cone $B$N7W;;$,I,MW(B.
% R^d $B$N(B cone $B$O(B cone.cinit.d $B$XF~$l$k(B.
cone.gw_ineq_projectedWtLpt getConeInfo /cone.cinit.d set
} ifelse
cone.cinit.d 1 get newCone /cone.startingCone set
(cone.startingCone is ) message
cone.startingCone message
] pop
popVariables
cone.startingCone
} def
%
% data/test9.sm1 $B$N(B test9 1-simplex X 2-simplex
%
% data/test10.sm1 1-simplex X 3-simplex
% data/test11.sm1 SST, p.59
%
% $B$$$h$$$h(B, cone enumeration $B$N%W%m%0%i%`=q$-3+;O(B
%
%<
% Usages: cone markBorder
% cone->facets[i] $B$,(B weight space $B$N(B border $B$K$"$k$H$-(B
% cone->flipped[i] = 2 $B$H$9$k(B.
% $B$3$l$r(B cone $B$N$9$Y$F$N(B facet $B$KBP$7$F7W;;(B.
%>
/markBorder {
/arg1 set
[/cone /facets_t /flipped_t /kk] pushVariables
[
/cone arg1 def
cone (facets) getNode 2 get /facets_t set
cone (flipped) getNode 2 get /flipped_t set
0 1 flipped_t length 1 sub {
/kk set
flipped_t kk get (0).. eq {
cone kk isOnWeightBorder {
% Border $B$N>e$K$"$k$N$G(B flip $B:Q$N%^!<%/$r$D$1$k(B.
flipped_t kk (2).. put
} { } ifelse
} { } ifelse
} for
] pop
popVariables
} def
%<
% Usages: ncone updateFan
% $B%0%m!<%P%kJQ?t(B cone.fan $B$r99?7$9$k(B.
%>
%
% updateFan $B$N(B debug $B$O(B data/test8 $B$G$H$j$"$($:$d$k(B.
% test8 /ncone set $B$r<B9T$7$F$+$i(B ncone updateFan
% global: cone.fan
/cone.fan [ ] def
% global: cone.incidence
/cone.incidence [ ] def
/updateFan {
/arg1 set
[/ncone /kk /cfacet /ii /jj /tcone /flipped_t] pushVariables
[
/ncone arg1 def
/cone.fan.n cone.fan length def
% 0. ncone $B$,(B cone.fan $B$K$9$G$K$"$l$P%(%i!<(B
0 1 cone.fan.n 1 sub {
/kk set
ncone cone.fan kk get isSameCone {
(Internal error updateFan: ncone is already in cone.fan) error
} { } ifelse
} for
% 1. ncone $B$NCf$N(B border $B>e$N(B facet $B$r$9$Y$F(B mark.
ncone markBorder
% 2. ncone /\ cone.fan[kk] $B$,$"$k$+D4$Y$k(B. $B$"$l$P(B Mark $B$9$k(B. incidence graph $B$K2C$($k(B
0 1 cone.fan.n 1 sub {
/kk set
ncone cone.fan kk get getCommonFacet /cfacet set
cfacet 0 get
{
% $B6&DL(B facet $B$,$"$k>l9g(B. [[cone$BHV9f(B face$BHV9f(B] [cone$BHV9f(B face$BHV9f(B]] $B$N7A<0$G3JG<(B.
/ii cfacet 1 get 0 get def
/jj cfacet 2 get 0 get def
cone.incidence [ [[cone.fan.n ii] [kk jj]] ] join /cone.incidence set
% flipped $B$r(B mark $B$9$k(B.
ncone ii markFlipped
cone.fan kk get /tcone set
tcone jj markFlipped
} { } ifelse
} for
% 3. ncone $B$r2C$($k(B.
cone.fan [ncone] join /cone.fan set
] pop
popVariables
} def
%<
% usages: getNextFlip [cone, k]
% cone.fan $B$r8!:w$7$F(B $B$^$@(B flip $B$7$F$J$$(B cone $B$H(B facet $B$NAH$rLa$9(B.
% $B$b$&$J$$$H$-$K$O(B null $B$rLa$9(B.
%>
/getNextFlip {
[/tcone /ans /ii ] pushVariables
[
/ans null def
0 1 cone.fan length 1 sub {
/ii set
cone.fan ii get /tcone set
tcone getNextFacet /ans set
ans tag 0 eq { } { exit } ifelse
} for
ans tag 0 eq { /arg1 null def }
{ /arg1 [tcone ans] def } ifelse
] pop
popVariables
arg1
} def
% global variable : cone.epsilon , cone.epsilon.limit
% flip $B$N;~$N(B epsilon
/cone.epsilon (1).. (10).. div def
/cone.epsilon.limit (1).. (100).. div def
%<
% Usages: result_getNextFlip getNextCone ncone
% flip $B$7$F?7$7$$(B ncone $B$rF@$k(B.
%>
/getNextCone {
/arg1 set
[/ncone /ccone /kk /w /next_weight_w_wv] pushVariables
[
/ccone arg1 def
/ncone null def
/kk ccone 1 get def
ccone 0 get /ccone set
{
ccone tag 0 eq { exit } { } ifelse
% ccone $B$N(B kk $BHVL\$N(B facet $B$K$D$$$F(B flip $B$9$k(B.
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.input next_weight_w_wv 1 get cone.gb /cone.cgb 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
% $B<!85$rD4$Y$k(B. $B$@$a$J$i(B 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 $B$,$J$1$l$P(B $B$d$O$j(B epsilon $B$r>.$5$/(B.
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
(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
%<
% Usages: set globals and getGrobnerFan
% cf. clearGlobals
% getStartingCone $B$9$k$H(B weightSpace $B$H$+$N7W;;$,$G$-$k(B. isOnWeightBorder $B$,(B
% $B7h$a$i$l$k(B.
%>
% $B$H$j$"$($:(B (data/test8.sm1) run $B$7$F$+$i(B getGrobnerFan
/getGrobnerFan {
getStartingCone /cone.ncone set
{
cone.ncone updateFan
( ) message
(----------------------------------------------------------) message
(getGrobnerFan #cone.fan=) messagen cone.fan length message
cone.ncone /cone.ccone set
getNextFlip /cone.nextflip set
cone.nextflip tag 0 eq { exit } { } ifelse
cone.nextflip getNextCone /cone.ncone set
} loop
(Construction is completed. See cone.fan and cone.incidence.) message
} def