File: [local] / OpenXM / src / hgm / mh / src / Testdata / test1.rr (download)
Revision 1.5, Tue Feb 9 05:00:31 2016 UTC (8 years, 4 months ago) by takayama
Branch: MAIN
CVS Tags: HEAD Changes since 1.4: +2 -2
lines
hgm_jack-n and hgm_jack-n-2f1 are unified to hgm_jack-n.
New input format (%!version2.0) is implemented.
Samples are Testdata/new-*-in.txt
|
/* $OpenXM: OpenXM/src/hgm/mh/src/Testdata/test1.rr,v 1.5 2016/02/09 05:00:31 takayama Exp $
*/
import("tk_jack.rr")$
import("ok_diff.rr")$
def test1() {
M=3;
Approx=9;
AA=[1/5,2/5]; BB=[3/7];
Z1=tk_jack.zonal([1],3);
printf("Approx=1: %a\n",F1=tk_jack.mh_t(AA,BB,M,1));
Z2=tk_jack.zonal([2],3);
Z11=tk_jack.zonal([1,1],3);
printf("Approx=2: %a\n",F2=tk_jack.mh_t(AA,BB,M,2));
F=tk_jack.mh_t([1/5,2/5],[3/7],M,Approx);
QK=[tk_jack.qk([1],AA,BB),tk_jack.qk([2],AA,BB),tk_jack.qk([1,1],AA,BB)];
X0g=0.166667;
Beta=[1,2,3];
Rule=[];
for (I=0; I<M; I++) {
Rule=cons([util_v(x,[I+1]),Beta[I]*X0g],Rule);
}
Rule=reverse(Rule);
Ans0=base_replace(F,Rule);
X0g=0.250360;
Beta=[1,2,3];
Rule=[];
for (I=0; I<M; I++) {
Rule=cons([util_v(x,[I+1]),Beta[I]*X0g],Rule);
}
Rule=reverse(Rule);
Ans1=base_replace(F,Rule);
}
// Equations for 2F1. cf. dlmf.nist.gov
def pp(Param,X,I) {
M=length(X);
A=Param[0]; B=Param[1]; C=Param[2];
Val=(C-(M-1)/2-(A+B+1-(M-1)/2)*X[I])/(X[I]*(1-X[I]));
return(Val);
}
def qq2(X,I,J) {
Val=(1/2)/(X[I]-X[J]);
return(Val);
}
def qq(X,I,J) {
Val=(1/2)*X[J]*(1-X[J]);
Val /= X[I]*(1-X[I])*(X[I]-X[J]);
return(Val);
}
def rr(Param,X,I) {
A=Param[0]; B=Param[1];
return(A*B/(X[I]*(1-X[I])));
}
def eq2f1_i(Param,X,I) {
M = length(X);
Dx = newvect(M);
for (J=0; J<M; J++) Dx[J] = util_v(dx,[J+1]);
Dx = vtol(Dx);
L = Dx[I]^2 + pp(Param,X,I)*Dx[I];
for (J=0; J<M; J++) {
if (J != I) {
L += qq2(X,I,J)*Dx[I];
}
}
for (J=0; J<M; J++) {
if (J != I) {
L += -qq(X,I,J)*Dx[J];
}
}
L += -rr(Param,X,I);
return(L);
}
def test2() {
if (type(getopt(approx))>0) {
Approx=getopt(approx);
}else Approx=4;
M=3;
X = newvect(M);
for (I=0; I<M; I++) X[I] = util_v(x,[I+1]);
X = vtol(X);
AA=[1/5,2/5]; BB=[3/7]; Param=append(AA,BB);
F=tk_jack.mh_t(AA,BB,M,Approx);
Rule=assoc(X,[t,2*t,3*t]);
R=newvect(M);
for (I=0; I<M; I++) {
R[I]=odiff_act(nm(eq2f1_i(Param,X,I)),F,X);
R[I]=red(base_replace(R[I],Rule));
}
printf("Approx=%a\n",Approx);
return(vtol(R));
}
// 2016.02.03
def gen3() {
Lam = [1,2,4,6,8,10];
M = length(Lam);
X = Lam[idiv(M,2)];
Beta = newvect(M);
for (I=0; I<M; I++) Beta[I] = X/(Lam[I]+X);
N1 = 5;
N2 = 10;
A=N1/2;
B=-N2/2+(M+1)/2;
C=(N1+M+1)/2;
X0=Beta[0]/10;
X1=1;
Param=[vtol(Beta),[A,B,C],X0,X1];
printf("%a\n",Param);
gen_input(Param);
}
// same parameter with Hashiguchi note.
def gen3b() {
Lam = [1,2,4];
M = length(Lam);
X = Lam[idiv(M,2)];
Beta = newvect(M);
for (I=0; I<M; I++) Beta[I] = X/(Lam[I]+X);
N1 = 5;
N2 = 10;
A=N1/2;
B=-N2/2+(M+1)/2;
C=(N1+M+1)/2;
X0=Beta[0]/10;
X1=1;
Param=[vtol(Beta),[A,B,C],X0,X1];
printf("[Beta,[A,B,C],X0,X1]=%a\n",Param);
gen_input(Param);
}
def gen_input(Param) {
if (type(getopt(fname)) > 0) {
Fname=rtostr(getopt(fname));
}else Fname="tmp-2f1-in.txt";
Beta=Param[0];
M = length(Beta);
A=Param[1][0];
B=Param[1][1];
C=Param[1][2];
X0=Param[2];
X1=Param[3];
Fp = open_file(Fname,"w");
fprintf(Fp,"%Mg=\n%a\n",M);
for (I=0; I<M; I++) {
fprintf(Fp,"%Beta[%a]=\n%a\n",I,deval(Beta[I]));
}
fprintf(Fp,"%Ng=\n1.0\n");
fprintf(Fp,"%X0g=\n%a\n",deval(X0));
fprintf(Fp,"%Iv's=\n*\n");
fprintf(Fp,"%Ef=\n*\n");
fprintf(Fp,"%Hg=\n0.001\n");
fprintf(Fp,"%Dp=\n1\n");
fprintf(Fp,"%Xng=\n%a\n",deval(X1));
fprintf(Fp,"%p_pFq=2, %a, %a\n",deval(A),deval(B));
fprintf(Fp,"%q_pFq=1, %a\n",deval(C));
close_file(Fp);
return(Fname);
}
// 2016.02.04. We do not use the Kummer relation cf. gen3b()
// It seems to be more stable to use this formula.
// By Hashiguchi formula (positive beta), degree seems to be 11
// By takemura formula (negative beta), degree must be 18
// gen3c and before are for old hgm_jack-n-2f1
def gen3c() {
Lam = [1,2,4];
M = length(Lam);
X = Lam[idiv(M,2)];
Beta = newvect(M);
// for (I=0; I<M; I++) Beta[I] = X/(Lam[I]+X);
for (I=0; I<M; I++) Beta[I] = -Lam[I];
N1 = 5;
N2 = 10;
A=(M+1)/2;
B=(N1+N2)/2;
C=(N1+M+1)/2;
// X0=Beta[0]/10;
X0=1/15;
X1=1;
Param=[vtol(Beta),[A,B,C],X0,X1];
printf("[Beta,[A,B,C],X0,X1]=%a\n",Param);
printf("../hgm_jack-n-2f1 --idata tmp-2f1-in.txt --degree ?? >tt-2016-02-04-3c.txt\n");
printf("../hgm_w-n-2f1 --idata tt-2016-02-04-3c.txt --verbose\n");
gen_input(Param);
}
// For new hgm_jack-n-2f1
def gen4() {
Lam = [1,2,4];
N1 = 5;
N2 = 10;
M = length(Lam);
X = Lam[idiv(M,2)];
Beta = newvect(M);
for (I=0; I<M; I++) Beta[I] = Lam[I];
if (N1 < M || N2 < M) debug("N1 or N2 is small.");
A=(M+1)/2;
B=(N1+N2)/2;
C=(N1+M+1)/2;
X0=Beta[0]/10;
X1=4;
Param=[vtol(Beta),[A,B,C],X0,X1];
printf("[Beta,[A,B,C],X0,X1]=%a\n",Param);
printf("../hgm_jack-n-2f1 --idata tmp-2f1-in.txt --degree ?? >tt-2016-02-04-4.txt\n");
printf("../hgm_w-n-2f1 --idata tt-2016-02-04-4.txt --verbose\n");
gen_input(Param);
}
def gen5() {
Lam = [1,2,3,4,5,6,7,8];
N1 = 10;
N2 = 15;
M = length(Lam);
X = Lam[idiv(M,2)];
Beta = newvect(M);
for (I=0; I<M; I++) Beta[I] = Lam[I];
if (N1 < M || N2 < M) debug("N1 or N2 is small.");
A=(M+1)/2;
B=(N1+N2)/2;
C=(N1+M+1)/2;
X0=1/3;
X1=80;
Param=[vtol(Beta),[A,B,C],X0,X1];
printf("[Beta,[A,B,C],X0,X1]=%a\n",Param);
printf("../hgm_jack-n-2f1 --idata tmp-2f1-in.txt --degree ?? >tt-2016-02-04-5.txt\n");
printf("../hgm_w-n-2f1 --idata tt-2016-02-04-5.txt [--gnuplotf test-g] --verbose\n");
printf("Change the plot range by hand in test-g-gp.txt\n");
gen_input(Param);
}
end$