/* $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$