2元分割表HGM関数

Risa/Asir 2元分割表HGM関数説明書

3.0 版

2019 年 6 月 7 日

by Y.Goto, Y.Tachibana, N.Takayama

Copyright © Risa/Asir committers 2004–2019. All rights reserved.

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]
1 2元分割表HGMの関数説明書について  
2 2元分割表HGMの関数   * modular計算 * Binary splitting
Index  

1 2元分割表HGMの関数説明書について

この説明書では HGM(holonomic gradient method) を用いた2元分割表の関数について説明する. ChangeLog の項目は www.openxm.org の cvsweb で ソースコードを読む時の助けになる情報が書かれている. このパッケージは下記のようにロードする. 

load("gtt_ekn3.rr");

gtt_ekn3.rr は gtt_ekn.rr を置き換える大きく改良されたパッケージである. 最新版の asir-contrib package を取得するには, 下記のように更新関数を呼び出す. 

import("names.rr");
asir_contrib_update(|update=1);

本文中で引用している文献を列挙する.

このマニュアルで説明する関数を用いたプログラム例は gtt_ekn3/test-t1.rr など.

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2 2元分割表HGMの関数

2.1.1 gtt_ekn3.gmvector  
2.1.2 gtt_ekn3.nc  
2.1.3 gtt_ekn3.lognc  
2.1.4 gtt_ekn3.expectation  
2.1.5 gtt_ekn3.setup  
• gtt_ekn3.upAlpha  
2.1.7 gtt_ekn3.cmle  
• gtt_ekn3.set_debug_level  
• gtt_ekn3.contiguity_mat_list_2  
• gtt_ekn3.show_path  
• gtt_ekn3.get_svalue  
• gtt_ekn3.assert1  
• gtt_ekn3.assert2  
• gtt_ekn3.assert3  
2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1  
[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1 超幾何関数E(k,n)

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.1 gtt_ekn3.gmvector

gtt_ekn3.gmvector(beta,p)

:: 周辺和 beta, セルの確率 p の二元分割表に付随する超幾何関数 E(k,n) の値およびその微分の値を戻す.

gtt_ekn3.ekn_cBasis_2(beta,p)

の別名である.

return

ベクトル, 超幾何関数の値とその微分. 詳しくは下記.

beta

行和, 列和のリスト. 成分はすべて正であること.

p

二元分割表のセルの確率のリスト

例: 次は2 x 2 分割表で行和が [5,1], 列和が [3,3], 各セルの確率が [[1/2,1/3],[1/7,1/5]] の場合の gmvector の値である. 

[3000] load("gtt_ekn3.rr");
[3001] gtt_ekn3.gmvector([[5,1],[3,3]],[[1/2,1/3],[1/7,1/5]])
[775/27783]
[200/9261]

例: N を2以上の自然数とする時, Gauss の超幾何関数(この場合は多項式となる) F(-36N,-11N,2N,(1-1/N)/56) の値は T3 に代入される ( [TGKT] ). 

N=2; 
T2=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],[[1,(1-1/N)/56],[1,1]])[0][0];
D=fac(36*N)*fac(11*N)*fac(2*N-1); 
T3=T2*D;

ちなみに同じ値を Mathematica に計算させるには 

n=2; Hypergeometric2F1[-36*n,-11*n,2*n,(1-1/n)/56]

例: interval option 

[4009] P=gtt_ekn3.prob1(5,5,100);
[[[100,200,300,400,500],[100,200,300,400,500]],[[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1/23,1/29,1/31,1/37],[1,1/41,1/43,1/47,1/53],[1,1,1,1,1]]]

[4010] util_timing(quote(gtt_ekn3.gmvector([[100,200,300,400,500],[100,200,300,400,500]], [[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1/23,1/29,1/31,1/37],[1,1/41,1/43,1/47,1/53],[1,1,1,1,1]])))[1]; 
[cpu,72.852,gc,0,memory,4462742364,real,72.856]

[4011] util_timing(quote(gtt_ekn3.gmvector([[100,200,300,400,500],[100,200,300,400,500]], [[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1/23,1/29,1/31,1/37],[1,1/41,1/43,1/47,1/53],[1,1,1,1,1]]|interval=100)))[1];
[cpu,67.484,gc,0,memory,3535280544,real,67.4844]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な 計算ができる. 守備範囲の異なるプログラム同士の比較, debug 用参考. 

[3080] import("tk_fd.rr");
[3081] A=tk_fd.marginal2abc([4,5],[2,4,3]);
[-4,[-4,-3],-1]  // 2変数 FD のパラメータ. a,[b1,b2],c
[3082] tk_fd.fd_hessian2(A[0],A[1],A[2],[1/2,1/3]);
Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
[4483/124416,[ 1961/15552 185/1728 ],
 [ 79/288 259/864 ]
 [ 259/864 47/288 ]]
// 戻値は [F=F_D, gradient(F), Hessian(F)]

// ekn_gt での例と同じパラメータ.
[3543] A=tk_fd.marginal2abc([5,1],[3,3]);
[-5,[-3],-1]
[3544] tk_fd.fd_hessian2(A[0],A[1],A[2],[(1/3)*(1/7)/((1/2)*(1/5))]);
Computing Dmat(ca) for parameters B=[-3],X=[ 10/21 ]
[775/27783,[ 20/147 ],[ 17/42 ]]

参考: 一般の A 分布の正規化定数についての Hessian の計算は実験的 package ot_hessian_ahg.rr で実装のテストがされている. (これはまだ未完成のテスト版なので出力形式等も将来的には変更される.) 

import("ot_hgm_ahg.rr");
import("ot_hessian_ahg.rr");
def  htest4() {
  extern C11_A;
  extern C11_Beta;
  Hess=newmat(7,7);
  A =C11_A;
  Beta0= [b0,b1,b2,b3];
  BaseIdx=[4,5,6];
  X=[x0,x1,x2,x3,x4,x5,x6];
  for (I=0; I<7; I++) for (J=0; J<7; J++) {
    Idx = [I,J];
    H=hessian_simplify(A,Beta0,X,BaseIdx,Idx);
    Hess[I][J]=H;
    printf("[I,J]=%a, Hessian_ij=%a\n",Idx,H);
  } 
  return(Hess);
}
[2917] C11_A;
[[0,0,0,1,1,1,1],[1,0,0,1,0,1,0],[0,1,1,0,1,0,1],[1,1,0,1,1,0,0]]
[2918] C11_Beta;
[166,36,290,214]
[2919] Ans=htest4$
[2920] Ans[0][0];
[[((b1-b0-1)*x4)/(x0^2),[4]],[((b1-b0-1)*x6)/(x0^2),[6]],
 [(b1^2+(-2*b0-1)*b1+b0^2+b0)/(x0^2),[]],[(x6)/(x0),[6,0]],[(x4)/(x0),[4,0]]]
参照

gtt_ekn3.setup @ref{gtt_ekn3.pfaffian_basis}

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.2 gtt_ekn3.nc

gtt_ekn3.nc(beta,p)

:: 周辺和 beta, セルの確率 p の二元分割表の条件付き確率の正規化定数 Z およびその微分の値を戻す.

return

ベクトル [Z,[[d_11 Z, d_12 Z, ...], ..., [d_m1 Z, d_m2 Z, ...., d_mn Z]]]

beta

行和, 列和のリスト. 成分はすべて正であること.

p

二元分割表のセルの確率のリスト

例: 2x3 分割表での Z とその微分の計算. 

[2237] gtt_ekn3.nc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
[4483/124416,[ 353/7776 1961/15552 185/1728 ]
[ 553/20736 1261/15552 1001/13824 ]]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な 計算ができる. 

[3076] import("tk_fd.rr");
[3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
[-4,[-4,-3],-1]
[3078] tk_fd.ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
[ 1 1 1 ]
Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
[4483/124416,[[353/7776,1961/15552,185/1728],
              [553/20736,1261/15552,1001/13824]]]
// 戻値は [Z, [[d_11 Z, d_12 Z, d_13 Z],
//             [d_21 Z, d_22 Z, d_23 Z]]] の値. 
//           ここで d_ij は i,j 成分についての微分を表す.
参照

gtt_ekn3.setup gtt_ekn3.lognc

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.3 gtt_ekn3.lognc

gtt_ekn3.lognc(beta,p)

:: 周辺和 beta, セルの確率 p の二元分割表の条件付き確率の正規化定数 Z の log の近似値およびその微分の近似値を戻す.

return

ベクトル [log(Z), [[d_11 log(Z), d_12 log(Z), ...], [d_21 log(Z),...], ... ]

beta

行和, 列和のリスト. 成分はすべて正であること.

p

二元分割表のセルの確率のリスト

例: 2 × 3 分割表での例. 第一成分のみ近似値. 

[2238] gtt_ekn3.lognc([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
[-3.32333832422461674630,[ 5648/4483 15688/4483 13320/4483 ]
[ 3318/4483 10088/4483 9009/4483 ]]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な 計算ができる. 

[3076] import("tk_fd.rr");
[3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
[-4,[-4,-3],-1]
[3078] tk_fd.log_ahmat_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
[ 1 1 1 ]
Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
[-3.32333832422461674639485797719209322217260539267246045320,
 [[1.2598706, 3.499442, 2.971224],
  [0.7401293, 2.250278, 2.009591]]]
// 戻値は [log(Z), 
//          [[d_11 log(Z), d_12 log(Z), d_13 log(Z)],
//           [d_21 log(Z), d_22 log(Z), d_23 log(Z)]]] 
// の近似値.
参照

gtt_ekn3.setup gtt_ekn3.nc

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.4 gtt_ekn3.expectation

gtt_ekn3.expectation(beta,p)

:: 周辺和 beta, セルの確率 p の二元分割表の期待値を計算する.

return

二元分割表の各セルの期待値のリスト.

beta

行和, 列和のリスト. 成分はすべて正であること.

p

二元分割表のセルの確率のリスト

2×2, 3×3 の分割表の期待値計算例. 

[2235] gtt_ekn3.expectation([[1,4],[2,3]],[[1,1/3],[1,1]]);
[ 2/3 1/3 ]
[ 4/3 8/3 ]
[2236] gtt_ekn3.expectation([[4,5],[2,4,3]],[[1,1/2,1/3],[1,1,1]]);
[ 5648/4483 7844/4483 4440/4483 ]
[ 3318/4483 10088/4483 9009/4483 ]

[2442] gtt_ekn3.expectation([[4,14,9],[11,6,10]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]]);
[ 207017568232262040/147000422096729819 163140751505489940/147000422096729819 
                                        217843368649167296/147000422096729819 ]
[ 1185482401011137878/147000422096729819 358095302885438604/147000422096729819 
                                         514428205457640984/147000422096729819 ]
[ 224504673820628091/147000422096729819 360766478189450370/147000422096729819 
                                        737732646860489910/147000422096729819 ]

参考: 2 x m 分割表(Lauricella FD)についてはパッケージ tk_fd でも下記のように同等な 計算ができる. 

[3076] import("tk_fd.rr");
[3077] A=tk_fd.marginal2abc([4,5],[2,4,3]);
[-4,[-4,-3],-1]
[3078] tk_fd.expectation_abc(A[0],A[1],A[2],[[1,1/2,1/3],[1,1,1]]);
RS=[ 4 5 ], CSnew=[ 2 4 3 ], Ynew=[ 1 1/2 1/3 ]
[ 1 1 1 ]
Computing Dmat(ca) for parameters B=[-4,-3],X=[ 1/2 1/3 ]
[[5648/4483,7844/4483,4440/4483],
 [3318/4483,10088/4483,9009/4483]]
// 各セルの期待値.

参考: 一般の A 分布の計算は ot_hgm_ahg.rr. まだ実験的なため, module 化されていない. ot_hgm_ahg.rr についての参考文献: K.Ohara, N.Takayama, Pfaffian Systems of A-Hypergeometric Systems II — Holonomic Gradient Method, arxiv:1505.02947 

[3237] import("ot_hgm_ahg.rr");
// 2 x 2 分割表.
[3238] hgm_ahg_expected_values_contiguity([[0,0,1,1],[1,0,1,0],[0,1,0,1]],
        [9,6,8],[1/2,1/3,1/5,1/7],[x1,x2,x3,x4]|geometric=1);
oohg_native=0, oohg_curl=1
[1376777025/625400597,1750225960/625400597,
 2375626557/625400597,3252978816/625400597]
// 2 x 2 分割表の期待値.

// 2 x 3 分割表.
[3238] hgm_ahg_expected_values_contiguity(
 [[0,0,0,1,1,1],[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1]],
 [5,2,4,3],[1,1/2,1/3,1,1,1],[x1,x2,x3,x4,x5,x6]|geometric=1);
[5648/4483,7844/4483,4440/4483,3318/4483,10088/4483,9009/4483]
// 2 x 3 分割表の期待値. 上と同じ問題.

3 x 3 分割表. 構造的0が一つ. 

/*
  dojo, p.221 のデータ.  成績3以下の生徒は集めてひとつに.
  2 1 1
  8 3 3
  0 2 6

  row sum: 4,14,8
  column sum: 10,6,10
  0 を一つ含むので, (3,6) 型の A から 7 列目を抜く.
*/

A=[[0,0,0,1,1,1, 0,0],
   [0,0,0,0,0,0, 1,1],
   [1,0,0,1,0,0, 0,0],
   [0,1,0,0,1,0, 1,0],
   [0,0,1,0,0,1, 0,1]];
B=[14,8,10,6,10];
hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1],
　　　　　　　　　　　　　　　　[x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);

// 答.
[14449864949304/9556267369631,
 10262588586540/9556267369631, 13512615942680/9556267369631,
 81112808747006/9556267369631,
 21816297744346/9556267369631, 30858636683482/9556267369631,
                              
 25258717886900/9556267369631,51191421070148/9556267369631]

3 x 3 分割表. 

/*
 上のデータで 0 を 1 に変更.
  2 1 1
  8 3 3
  1 2 6

  row sum: 4,14,9
  column sum: 11,6,10
*/
A=[[0,0,0,1,1,1,0,0,0],
   [0,0,0,0,0,0,1,1,1],
   [1,0,0,1,0,0,1,0,0],
   [0,1,0,0,1,0,0,1,0],
   [0,0,1,0,0,1,0,0,1]];
B=[14,9,11,6,10];
hgm_ahg_expected_values_contiguity(A,B,[1,1/2,1/3,1,1/5,1/7,1,1,1],
                              [x1,x2,x3,x4,x5,x6,x7,x8]|geometric=1);

// 期待値, 答.   x9 を指定していないので, 9番目の期待値は出力してない.
[207017568232262040/147000422096729819,
 163140751505489940/147000422096729819,217843368649167296/147000422096729819,
 1185482401011137878/147000422096729819,
 358095302885438604/147000422096729819,514428205457640984/147000422096729819,
 224504673820628091/147000422096729819,360766478189450370/147000422096729819]

// Z やその微分の計算は hgm_ahg_contiguity 関数がおこなうが, これの簡易インターフェースは
// まだ書いてない.
参照

gtt_ekn3.setup gtt_ekn3.nc

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.5 gtt_ekn3.setup

gtt_ekn3.setup()

:: 分散計算用の環境設定をおこなう. 現在の環境を報告する.

return

例: 素数のリストを生成してファイル p.txt へ書き出す. 

gtt_ekn3.setup(|nps=2,nprm=20,minp=10^10,fgp="p.txt")$

例: chinese remainder theorem (crt) を使って gmvector を計算. 

[2867] gtt_ekn3.setup(|nprm=20,minp=10^20);
[2868] N=2; T2=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                                [[1,(1-1/N)/56],[1,1]] | crt=1)$
参照

gtt_ekn3.nc gtt_ekn3.gmvector

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.6 gtt_ekn3.upAlpha, gtt_ekn3.downAlpha

gtt_ekn3.upAlpha(i,k,n)
gtt_ekn3.downAlpha(i,k,n)

::

i a_i を a_i+1 (a_i を a_i-1) と変化させる contiguity relation.
k E(k+1,n+k+2)型の超幾何関数の k. 分割表では (k+1)×(n+1).
n E(k+1,n+k+2)型の超幾何関数の n. 分割表では (k+1)×(n+1).
return contiguity relation の pfaffian_basis についての行列表現を戻す. [GM2016] の Cor 6.3.

例: 以下の例は 2×2分割表(E(2,4)), 2×3分割表(E(2,5))の場合である. [2225] までは出力を略している. 

[2221] gtt_ekn3.marginaltoAlpha([[1,4],[2,3]]);
[[a_0,-4],[a_1,-1],[a_2,3],[a_3,2]]
[2222] gtt_ekn3.upAlpha(1,1,1);  // E(2,4) の a_1 方向の 
                                //     contiguity を表現する行列
[2223] gtt_ekn3.upAlpha(2,1,1);  // E(2,4) の a_2 方向
[2224] gtt_ekn3.upAlpha(3,1,1);  // E(2,4) の a_3 方向
[2225] function f(x_1_1);
[2232] gtt_ekn3.pfaffian_basis(f(x_1_1),1,1);
[ f(x_1_1) ]
[ (f{1(x_1_1)*x_1_1)/(a_2) ]
[2233] function f(x_1_1,x_1_2);
f() redefined.
[2234] gtt_ekn3.pfaffian_basis(f(x_1_1,x_1_2),1,2); // E(2,5), 2*3 分割表
[ f(x_1_1,x_1_2) ]
[ (f{1,0(x_1_1,x_1_2)*x_1_1)/(a_2) ]
[ (f{0,1(x_1_1,x_1_2)*x_1_2)/(a_3) ]

[2235]   RuleA=[[a_2,1/3],[a_3,1/2]]$ RuleX=[[x_1_1,1/5]]$
  base_replace(gtt_ekn3.upAlpha(1,1,1),append(RuleA,RuleX))
 -gtt_ekn3.upAlpha(1,1,1 | arule=RuleA, xrule=RuleX);

[ 0 0 ]
[ 0 0 ]
参照

gtt_ekn3.nc gtt_ekn3.gmvector

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.7 gtt_ekn3.cmle

gtt_ekn3.cmle(u) u を観測データとするとき, P(U=u | row sum, column sum = these of U) を最大化する, 各セルの確率の近似値を求める.

::

u 観測データ(分割表)
return セルの確率(分割表形式)

例: 2 x 4 分割表. 

U=[[1,1,2,3],[1,3,1,1]];
gtt_ekn3.cmle(U);
 [[ 1 1 2 3 ]
  [ 1 3 1 1 ],[[7,6],[2,4,3,4]],   // Data, row sum, column sum
 [ 1 67147/183792 120403/64148 48801/17869 ]  // probability obtained.
 [ 1 1 1 1 ]]

例: 上の例は次の関数に. 

gtt_ekn3.cmle_test3();
参照

gtt_ekn3.expectation

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1

gtt_ekn3.set_debug_level(m) debug メッセージのレベルを設定.
gtt_ekn3.contiguity_mat_list_2 使用する contiguity を構成.
gtt_ekn3.show_path() どのように contiguity を適用したかの情報.
gtt_ekn3.get_svalue() static 変数の値を得る.
gtt_ekn3.assert1(N) 2x2 分割表で動作チェック.
gtt_ekn3.assert2(N) 3x3 分割表で動作チェック.
gtt_ekn3.assert3(R1, R2, Size) R1 x R2 分割表で並列動作の動作チェック.
gtt_ekn3.prob1(R1,R2,Size) R1 x R2 分割表用のテストデータを作る.

::

m レベル.

例. 

[2846] gtt_ekn3.set_debug_level(0x4);
[2847] N=2; T2=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                                [[1,(1-1/N)/56],[1,1]])$
[2848] level&0x4: g_mat_fac_test([ 113/112 ]
[ 1/112 ],[ (t+225/112)/(t^2+4*t+4) (111/112*t+111/112)/(t^2+4*t+4) ]
[ (1/112)/(t^2+4*t+4) (111/112*t+111/112)/(t^2+4*t+4) ],0,20,1,t)
Note: we do not use g_mat_fac_itor. Call gtt_ekn3.setup(); to use the crt option.
level&0x4: g_mat_fac_test([ 67/62944040755546030080000 ]
[ 1/125888081511092060160000 ],[ (t+24)/(t^2+25*t+46) (2442)/(t^2+25*t+46) ]
[ (1)/(t^2+25*t+46) (-111*t-111)/(t^2+25*t+46) ],0,73,1,t)
level&0x4: g_mat_fac_test ------  snip

例. 

[2659] gtt_ekn3.nc([[4,5,6],[2,4,9]],[[1,1/2,1/3],[1,1/5,1/7],[1,1,1]])$
[2660] L=matrix_transpose(gtt_ekn3.show_path())$
[2661] L[2];
[2 1]

[2 1] の index をもつパラメーター alpha の方向の contigity を求めそれを掛けて 計算したことがわかる. L[0] は用いた contiguity の行列. L[1] はcontiguity を適用する step 数.

例. 値を計算せずに path のみ求めたい場合. 

A=gtt_ekn3.marginaltoAlpha_list([[400,410,1011],[910,411,500]])$
[2666] gtt_ekn3.contiguity_mat_list_2(A,2,2)$
[2667] L=matrix_transpose(gtt_ekn3.show_path())$
[2668] L[2];
[ 2 1 5 4 3 ]
[2669] gtt_ekn3.contiguity_mat_list_3(A,2,2)$ // new alg in [TGKT]
[2670] L=matrix_transpose(gtt_ekn3.show_path())$
[2671] L[2];
[2 1]  // shorter

例. 値を計算せずに path のみ求めたい場合. gtt_ekn3 による新しいアルゴリズムによる path の表示. 

A=gtt_ekn3.marginaltoAlpha_list([[10,20],[15,15]])$
[2666] gtt_ekn3.contiguity_mat_list_3(A,1,1 | xrule=[[x_1_1,1/2]])$
[t,[[ (-t-43/2)/(t-2) (-15/2)/(t-2) ]
[ 1/2 -1/2 ],-9]]

例. 0 が戻れば g_mat_fac_plain と指定した計算方法の結果が一致したことがわかる. option を書かないと g_mat_fac_int との比較となる. 

[8859] gtt_ekn3.assert2(1);
Marginal=[[130,170,353],[90,119,444]]
P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]
Try g_mat_fac_test_int: Note: we do not use g_mat_fac_itor. Call gtt_ekn3.setup(); to use the crt option.
Timing (int) =0.413916 (CPU) + 0.590723 (GC) = 1.00464 (total), real time=0.990672

Try g_mat_fac_test_plain: Note: we do not use g_mat_fac_itor. Call gtt_ekn3.setup(); to use the crt option.
Timing (rational) =4.51349 (CPU) + 6.32174 (GC) = 10.8352 (total)
diff of both method = 
[ 0 0 0 ]
[ 0 0 0 ]
[ 0 0 0 ]
[8860] 

[8863] gtt_ekn3.setup(|nprm=100,minp=10^50);
Number of processes = 1.
Number of primes = 100.
Min of plist = 100000000000000000000000000000000000000000000000151.
0
[8864] gtt_ekn3.assert2(1 | crt=1);
Marginal=[[130,170,353],[90,119,444]]
P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]
Try [[crt,1]]
----  snip

なお二番目の例の timing (total) [例では省略] は mod 計算を subprocess がやっているので正しい値ではない. real time が計算時間の目安になる.

例. 

3x5 分割表. 周辺和は 10 に比例する一定の数(factor option も関係. ソースを参照). 
cell 確率は1/素数で生成される.
[9054] L=gtt_ekn3.prob1(3,5,10 | factor=1, factor_row=3); 
[[[10,20,420],[30,60,90,120,150]],[[1,1/2,1/3,1/5,1/7],[1,1/11,1/13,1/17,1/19],[1,1,1,1,1]]]
[9055] number_eval(gtt_ekn3.expectation(L[0],L[1]));
[ 1.65224223218613 ... snip ]

例: 

[5779] import("gtt_ekn3.rr"); load("gtt_ekn3/ekn_eval-timing.rr");
[5780] gtt_ekn3.assert3(5,5,100 | nps=32, interval=100);
 -- snip
Parallel method: Number of process=32, File name tmp-gtt_ekn3/p300.txt is written.
Number of processes = 32.
  -- snip
initialPoly of path=3: [ 2.184 0 124341044 2.1831 ] [CPU(s),0,*,real(s)]
contiguity_mat_list_3 of path=3: [ 0.04 0 630644 9.6774 ] [CPU(s),0,*,real(s)]
Note: interval option will lead faster evaluation. We do not use g_mat_fac_itor (crt). Call gtt_ekn3.setup(); to use the crt option.
g_mat_fac of path=3: [ 21.644 0 1863290168 21.6457 ] [CPU(s),0,*,real(s)]
Done. Saved in 2.ab
Diff (should be 0)=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,..., 0,0,0]
参照

gtt_ekn3.nc

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

3 modular計算

3.1.1 gtt_ekn3.chinese_itor  
[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

3.1 中国剰余定理とitor

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

3.1.1 gtt_ekn3.chinese_itor

gtt_ekn3.chinese_itor(data,idlist)

:: mod p で計算した結果(ベクトル)から chinese remainder theorem, itor(integer to rational) で有理数ベクトルを得る.

return [val, n] ここで val は答え. また, n = n1*n2*...
data [[val1,n1],[val2,n2], ...], ここで val mod n1 = val1, val mod n2 = val2,...
idlist chinese, itor を実行するサーバIDのリスト.

例: [3!, 5^3*3!]=[6,750] が戻り値. 6 mod 109 =6, 750 mod 109=96 が最初の引数の [[6,96],109]. 以下同様. 

gtt_ekn3.setup(|nps=2,nprm=3,minp=101,fgp="p_small.txt");
SS=gtt_ekn3.get_svalue();
SS[0];
  [103,107,109]   // list of primes
SS[1];
  [0,2]           // list of server ID's
gtt_ekn3.chinese_itor([[[ 6,96 ],109],[[ 6,29 ],103],[[ 6,1 ],107]],SS[1]);
  [[ 6 750 ],1201289]

// 引数はスカラーでもよい.
gtt_ekn3.chinese_itor([[96,109],[29,103]],SS[1]);
  [[ 750 ],11227]

例: gtt_ekn3/childprocess.rr (server で実行される) の関数 chinese (chinese remainder theorem) と euclid. 

load("gtt_ekn3/childprocess.rr");
chinese([newvect(2,[6,29]),103],[newvect(2,[6,750]),107*109]);
  // mod 103 で [6,29], mod (107*109) で [6,750] となる数を mod 103*(107*109)
  // で求めると,
  [[ 6 750 ],1201289]
euclid(3,103);  // mod 103 での 3 の逆数. つまり 1/3
  -34
3*(-34) % 103; // 確かに逆数.
   1

例: gtt_ekn3/childprocess.rr (server で実行される) の関数 itor (integer to rational) の例. itor(Y,Q,Q2,Idx) では Y < Q2 なら Y がそのまま戻る. Idx は 内部用の index で好きな数でよい. 戻り値の第２成分となる. 

load("gtt_ekn3/childprocess.rr");
for (I=1;I<11; I++) print([I,itor(I,11,3,0)]);
[1,[1,0]]
[2,[2,0]]
[3,[-2/3,0]] //euclid(3,11); ->4,  4*(-2)%11 -> 3 なので確かに -2/3 は元の数の候補
[4,[failure,0]]
[5,[-1/2,0]]
[6,[1/2,0]]
[7,[-1/3,0]]
[8,[failure,0]]
[9,[-2,0]]
[10,[-1,0]]
参照

gtt_ekn3.setup

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

4 Binary splitting

• gtt_ekn3.init_dm_bsplit  
4.1.1 gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit  
• gtt_ekn3.init_bsplit  
[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

4.1 matrix factorial

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

4.1.1 gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit

gtt_ekn3.init_bsplit(|minsize=16,levelmax=1);

:: binary split の実行のためのパラメータを設定する.

gtt_ekn3.init_dm_bsplit(|bsplit_x=0, bsplit_reduce=0)

:: binary split の分散実行のためのパラメータを設定する.

gtt_ekn3.setup_dm_bsplit(C)

:: binary split の分散実行のために C 個のプロセスを立ち上げる.

C はlevelmax-1 に設定する. 特に levalmax=1 のときは分散計算を行わない.
bsplit_x=1 のとき, debug 用に各プロセスを xterm で表示.

例: bs=1 と無い場合の比較. 

[4618] cputime(1)$
[4619] gtt_ekn3.expectation(Marginal=[[1950,2550,5295],[1350,1785,6660]],
                          P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]]|bs=1)$
4.912sec(4.914sec)
[4621] V2=gtt_ekn3.expectation(Marginal=[[1950,2550,5295],[1350,1785,6660]],
                          P=[[17/100,1,10],[7/50,1,33/10],[1,1,1]])$
6.752sec(6.756sec)

例: 分散計算する場合. 分散計算はかえって遅くなる場合が多いので注意. 下記の例での bsplit_x=1 option は debug windows を開くのでさらに遅くなる. gtt_ekn3.test_bs_dist(); でもテストできる. 

[3669] C=4$ gtt_ekn3.init_bsplit(|minsize=16,levelmax=C+1)$ gtt_ekn3.init_dm_bsplit(|bsplit_x=1)$
[3670] [3671] [3672] gtt_ekn3.setup_dm_bsplit(C);
[0,0]
[3673] gtt_ekn3.assert2(10|bs=1)$
参照

gtt_ekn3.gmvector gtt_ekn3.expectation gtt_ekn3.assert1 gtt_ekn3.assert2

ChangeLog

[ << ] [ < ] [上] [ > ] [ >> ]         [冒頭] [目次] [見出し] [ ? ]

Index

移動:   G  
見出し一覧  項
G
gtt_ekn3.assert1 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.assert2 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.assert3 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.chinese_itor 中国剰余定理とitor 3.1.1 gtt_ekn3.chinese_itor
gtt_ekn3.cmle 2.1.7 gtt_ekn3.cmle
gtt_ekn3.contiguity_mat_list_2 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.downAlpha 2.1.6 gtt_ekn3.upAlpha, gtt_ekn3.downAlpha
gtt_ekn3.expectation 2.1.4 gtt_ekn3.expectation
gtt_ekn3.get_svalue 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.gmvector 2.1.1 gtt_ekn3.gmvector
gtt_ekn3.init_bsplit matrix factorial 4.1.1 gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit
gtt_ekn3.init_dm_bsplit matrix factorial 4.1.1 gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit
gtt_ekn3.lognc 2.1.3 gtt_ekn3.lognc
gtt_ekn3.nc 2.1.2 gtt_ekn3.nc
gtt_ekn3.prob1 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.setup 2.1.5 gtt_ekn3.setup
gtt_ekn3.setup_dm_bsplit matrix factorial 4.1.1 gtt_ekn3.init_bsplit, gtt_ekn3.init_dm_bsplit, gtt_ekn3.setup_dm_bsplit
gtt_ekn3.set_debug_level 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.show_path 2.1.8 gtt_ekn3.set_debug_level, gtt_ekn3.show_path, gtt_ekn3.get_svalue, gtt_ekn3.assert1, gtt_ekn3.assert2, gtt_ekn3.assert3, gtt_ekn3.prob1
gtt_ekn3.upAlpha 2.1.6 gtt_ekn3.upAlpha, gtt_ekn3.downAlpha
移動:   G  
[冒頭] [目次] [見出し] [ ? ]

目次

[冒頭] [目次] [見出し] [ ? ]

簡略化した目次

[冒頭] [目次] [見出し] [ ? ]

この文書について

この文書は4月 20, 2024にtexi2html 5.0を用いて生成されました。

ナビゲーションパネル中のボタンには以下の意味があります。

ボタン 名称 移動先 1.2.3項からの移動先
[ << ] FastBack Beginning of this chapter or previous chapter 1
[ < ] Back Previous section in reading order 1.2.2
[上] Up Up section 1.2
[ > ] Forward Next section in reading order 1.2.4
[ >> ] FastForward Next chapter 2
[冒頭] 冒頭 Cover (top) of document  
[目次] 目次 Table of contents  
[見出し] 見出し 見出し  
[ ? ] About About (help)  

例では、以下に示す構造を持つ文書の1.2.3項を現在位置に仮定しています。

この文書は4月 20, 2024にtexi2html 5.0を用いて生成されました。