Return to gtt_ekn3.rr CVS log

Up to [local] / OpenXM / src / asir-contrib / packages / src

gtt_ekn3.prob5x5(N) returns the benchmark problem 3 of [TGKT].

/* $OpenXM: OpenXM/src/asir-contrib/packages/src/gtt_ekn3.rr,v 1.12 2019/06/18 01:32:32 takayama Exp $ */
import("names.rr")$
load("tk_nm_dn.rr")$
if ((type(version(2018)) == 4) && (length(version(2018)) < 3)) {
 igcdcntl(3)$ printf("igcdcntl is set to 3.\n")$
 printf("Warning: it is strongly recommended to run this package on asir-2018.\n")$
}else {;}

#define USE_MODULE

#ifdef USE_MODULE
module gtt_ekn3$
/* prototype declaration of all functions */
localf abs_max$ 
localf alpha$ 
localf alphaProd$ 
localf alphaRule_num$ 
localf alphaRule$ 
localf antiDiag$ 
localf append_comp_v$ 
localf assert1$ 
localf assert2$ 
localf assert3$ 
localf avoid_var$ 
localf basistodS_g_1$ 
localf basistodS$ 
localf basistoFdiff$ 
localf basistoHessian$ 
localf basistoS1_g_1$ 
localf basistoS1$ 
localf calJ$ 
localf calJ0_slow$ 
localf calJ0$ 
localf calJqp$ 
localf cBasistoE_0$ 
localf cBasistoE_g_1$ 
localf cBasistoE$ 
localf cBasistoEdiff_0$ 
localf cBasistoEdiff_g_1$ 
localf cBasistoEdiff$ 
localf check_failure$ 
localf check_relation$ 
localf chinese_itor$ 
localf cInt_g_1$ 
localf cInt$ 
localf cInts$ 
localf cInts2$ 
localf cmle_1$ 
localf cmle_2$ 
localf cmle_3_v1$ 
localf cmle_3$ 
localf cmle_options$ 
localf cmle_test1$ 
localf cmle_test1c$ 
localf cmle_test2$ 
localf cmle_test2b$ 
localf cmle_test2c$ 
localf cmle_test3$ 
localf cmle_test4$ 
localf cmle_test5_ans$ 
localf cmle_test5_female$ 
localf cmle_test5_male$ 
localf cmle_test5$ 
localf cmle_test5b$ 
localf cmle$ 
localf coef_upAlpha$ 
localf coefPsi$ 
localf contiguity_mat_list_1$ 
localf contiguity_mat_list_2$ 
localf contiguity_mat_list_3$ 
localf contiguity_mat_list_3_dist;
localf contiguity_mat_list_g_1$ 
localf contiguity_mat_list_g_2$ 
localf dataset$ 
localf denormalize$ 
localf denormalize2$ 
localf dg_mat_fac_mm$ 
localf dg_mat_fac_mm0$ 
localf dlogxJ_g_1$ 
localf dlogxJ$ 
localf dnlcm$ 
localf do_itor$ 
localf downAlpha_g_1$ 
localf downAlpha_g_2$ 
localf downAlpha$ 
localf downAlpha2$ 
localf downAlpha3$ 
localf downAlpha3_ax;
localf dxij$ 
localf e36_deg_contig_test_2$ 
localf e36_deg_contig_test_3$ 
localf e36_deg_contig_test_k$ 
localf e36_deg_contig_test_k2$ 
localf e36_deg_contig_test$ 
localf e36_deg_test_invintMatrix$ 
localf e36_deg_test$ 
localf e36_deg_timing_test$ 
localf e48_deg_contig_test$ 
localf ekn_cBasis_1$ 
localf ekn_cBasis_2$ 
localf ekn_cBasis_3$ 
localf ekn_cBasis_g_1$ 
localf ekn_restr_1$ 
localf ekn_restr_factor_1$ 
localf ell$ 
localf ellAll$ 
localf ellRule$ 
localf eulerxij$ 
localf eulerZ$ 
localf expectation$ 
localf factorial_prod$ 
localf factorial$ 
localf find_zero$ 
localf first_diff$ 
localf g_mat_fac_int$ 
localf g_mat_fac_itor$ 
localf g_mat_fac_test_bs$ 
localf g_mat_fac_test_int$ 
localf g_mat_fac_test_int2$ 
localf g_mat_fac_test_plain$ 
localf g_mat_fac_test$ 
localf generate_p_set_$ 
localf generate_p_set$ 
localf get_dm_bsplit_global$ 
localf get_ekn_IDL$
localf get_M_ans$
localf get_mfac_bs_global$ 
localf get_mfac_bs_global$ 
localf get_svalue$ 
localf get_T_list$ 
localf gm_test1$ 
localf gm_test1b$ 
localf gm_test2$ 
localf gm_test3$ 
localf gm_test3a$ 
localf gm_test3b$ 
localf gm_test3c$ 
localf gm_test4$ 
localf gmvector_0$ 
localf gmvector_g_1$ 
localf gmvector$ 
localf go_next$ 
localf has_negative_elem$ 
localf hpq$ 
localf ijphiJ$ 
localf init_bsplit$ 
localf init_dm_bsplit$ 
localf init_myrand$ 
localf initial_row$ 
localf initialAlpha_1$ 
localf initialAlpha_3$ 
localf initialAlpha$ 
localf initialBeta_1$ 
localf initialBeta_3$ 
localf initialBeta$ 
localf initialPoly_3$ 
localf initialPoly_table_3$ 
localf initialPoly_table$ 
localf initialPoly$ 
localf intMatrix_g_1$ 
localf intMatrix$ 
localf intMatrix2$ 
localf intMatrixH$ 
localf invintMatrix_g_1$ 
localf invintMatrix_g_l00l$ 
localf invintMatrix_k$ 
localf invintMatrix_k2$ 
localf invintMatrix_n$ 
localf invrepMatrix_g_1$ 
localf invrepMatrix$ 
localf invrepMatrix2$ 
localf is_no_zero_cell$ 
localf isConst$ 
localf jinv_orig$ 
localf jinv_v2$ 
localf jinv$ 
localf jj$ 
localf load_prime$ 
localf log_cont$ 
localf lognc$ 
localf ltop$ 
localf ltopoly$ 
localf lvec$ 
localf m_alpha$ 
localf marginaltoAlpha_list$ 
localf marginaltoAlpha$ 
localf mat_red$ 
localf max$ 
localf mBasis_replace_g_1$ 
localf mBasis_replace$ 
localf mBasis$ 
localf mBasisList$ 
localf mlcm$ 
localf mle_init$ 
localf modp_itor$ 
localf move$ 
localf moveAlpha_n$ 
localf moveAlpha$ 
localf mul_list$ 
localf myabs$ 
localf myrand$ 
localf mytiming$
localf nablaij$ 
localf nc$ 
localf next_2$ 
localf next$ 
localf nextData_itor$ 
localf nextP_itor$ 
localf nm_dn_assert$ 
localf nm_dn$ 
localf normalize$ 
localf number_nm_dn_0$ 
localf number_nm_dn_list$ 
localf number_nm_dn$ 
localf oddsvect$ 
localf pf_basis_zero$ 
localf pfaffian_basis_g_1$ 
localf pfaffian_basis$ 
localf pfaffianCoef_g_1$ 
localf pfaffianCoef$ 
localf pfaffianMJ$ 
localf pfaffianPsi_2$ 
localf pfaffianPsi_time$ 
localf pfaffianPsi_time2$ 
localf pfaffianPsi$ 
localf phiJ$ 
localf phitov$ 
localf pm$ 
localf pochhammer$ 
localf poly_nm_dn_list$ 
localf poly_nm_dn$ 
localf preparation_F$ 
localf preparation_Mat$ 
localf preparation$ 
localf prime_prob$ 
localf prime_XRule$ 
localf print_restr_factor_1$ 
localf prob1$ 
localf prob5x5$ 
localf process$ 
localf rdeg$ 
localf repMatrix_g_1$ 
localf repMatrix$ 
localf repMatrix2$ 
localf report$ 
localf rm_list$ 
localf s36$ 
localf s36$ 
localf save_prime$ 
localf set_ARule$ 
localf set_assert$ 
localf set_debug_level$ 
localf set_XRule$ 
localf setup_dm_bsplit$ 
localf setup_x$ 
localf setup$ 
localf shift_m$ 
localf shift$ 
localf show_path$ 
localf shutdown$ 
localf sortSgn$ 
localf submatrix$ 
localf sum_time$ 
localf swap_col$ 
localf swap_row$ 
localf test_bs_dist$ 
localf test5x5$
localf test_nxn$ 
localf test111$ 
localf test222$ 
localf ti$ 
localf to_vect$ 
localf tt$ 
localf tvec$ 
localf txmat$ 
localf upAlpha_g_1$ 
localf upAlpha_g_2$ 
localf upAlpha$ 
localf upAlpha2$ 
localf upAlpha3$ 
localf ux$ 
localf vanish_dvar$ 
localf vanish_var_ind$ 
localf vanish_var_list$ 
localf vanish_var$ 
localf vanish_vartoform$ 
localf xij_rule$ 
localf xij$ 
localf xJ_rule$ 
localf xJ_rule2$ 
localf xJ$ 
localf xmat$ 
localf xRule_num$

localf get_time_initial_poly$ // define in ekn_eval-timing.rr
localf get_time_contiguity$ 
localf get_time_g_mat_fac$ 
localf crt_debug$ // defined in g_mat_fac-debug.rr
localf crt_time_reset$
localf crt_time_get$
endmodule$
#else
extern Ekn_plist$
extern Ekn_IDL$
#endif
Gtt_ekn3="gtt_ekn3/"$
load(Gtt_ekn3+"g_mat_fac.rr")$
load(Gtt_ekn3+"dm_bsplit.rr")$
load(Gtt_ekn3+"mfac_bs.rr")$
load(Gtt_ekn3+"ekn_pfaffian_8.rr")$
load(Gtt_ekn3+"ekn_eval.rr")$
load(Gtt_ekn3+"mle.rr")$
load(Gtt_ekn3+"ekn_hgm_0.rr")$
load(Gtt_ekn3+"ekn_restriction_2.rr")$

#ifdef USE_MODULE
module gtt_ekn3;
#endif
def get_svalue() {  // get values of static or extern variables.
  return([Ekn_plist,Ekn_IDL,Ekn_debug,Ekn_mesg,XRule,ARule,Verbose,Ekn_Rq]);
}

// Gauss Manin vector
def gmvector(Beta,Prob) {
  if ((base_memberq(0,Beta[0])) || (base_memberq(0,Beta[1]))){
    error("gmvector: Marginal sums must be positive integers.");}
  if (base_sum(Beta[0],0,0,length(Beta[0])-1) != base_sum(Beta[1],0,0,length(Beta[1])-1)){
    error("gmvector: The sums of Beta[0] and Beta[1] must be equal.");}
  if (type(getopt(path)) < 0) Path=3;
  else if (getopt(path)<3) Path=2;
  else Path=3;
  if (Path == 2) return(ekn_cBasis_2(Beta,Prob|option_list=getopt()));
  else return(ekn_cBasis_3(Beta,Prob|option_list=getopt()));
}
def gmvector_g_1(Beta,Prob) {
  if ((base_memberq(0,Beta[0])) || (base_memberq(0,Beta[1]))){
    error("gmvector: Marginal sums must be positive integers.");}
  if (base_sum(Beta[0],0,0,length(Beta[0])-1) != base_sum(Beta[1],0,0,length(Beta[1])-1)){
    error("gmvector: The sums of Beta[0] and Beta[1] must be equal.");}
  return(ekn_cBasis_g_1(Beta,Prob|option_list=getopt()));
}

// Normalizing constant and its derivatives
def nc(Beta,Prob) {
  R1=length(Beta[0]); 
  R2=length(Beta[1]);
  GM=gmvector(Beta,Prob|option_list=getopt());
  NC=ltopoly(GM[0][0],Beta,Prob);
  Fdiff=base_replace(basistoFdiff(GM,R1-1,R2-1|xrule=xRule_num(Prob,R1-1,R2-1)),marginaltoAlpha(Beta));
  NCdiff=newmat(R1,R2);
  for (I=0; I<R1; I++){
    for (J=0; J<R2; J++){
      NCdiff[I][J]=red(eulerZ([I,J],Beta,Prob,GM[0][0],Fdiff))/Prob[I][J];
    }
  }  
  return([NC,NCdiff]);
}

// log of normalizing constant and its derivatives
def lognc(Beta,Prob) {
  R1=length(Beta[0]); 
  R2=length(Beta[1]);
  GM=gmvector(Beta,Prob | option_list=getopt());
  NC=ltopoly(GM[0][0],Beta,Prob);
  Fdiff=base_replace(basistoFdiff(GM,R1-1,R2-1|xrule=xRule_num(Prob,R1-1,R2-1)),marginaltoAlpha(Beta));
  LNCdiff=newmat(R1,R2);
  for (I=0; I<R1; I++){
    for (J=0; J<R2; J++){
      LNCdiff[I][J]=red(eulerZ([I,J],Beta,Prob,GM[0][0],Fdiff))/(Prob[I][J]*NC);
    }
  }  
  return([eval(log(NC)),LNCdiff]);
}

// Expectation for the index 
// expectation(Beta,Prob | index=[[0,0],[1,0]])
def expectation(Beta,Prob) {
  R1 = length(Beta[0]);
  R2 = length(Beta[1]);
  GM=gmvector(Beta,Prob | option_list=getopt());
  Expect=cBasistoE(GM,Beta,Prob);
  if (type(getopt(index)) >= 0) {
    IndexList = getopt(index);
    Ex=newmat(R1,R2);
    for (I=0; I<R1; I++){
      for (J=0; J<R2; J++){
	if (IndexList[I][J]==0){Ex[I][J]='*';}
	else{Ex[I][J]=Expect[I][J];}
      }
    }  
  }else{
    // Default IndexList
    IndexList=newmat(R1,R2);
    for (I=0; I<R1; I++) {
      for (J=0; J<R2; J++) {
	// IndexList[I*R2+J] = [I,J]; // ???
	IndexList[I][J] = 1;
      }
    }
    IndexList=matrix_matrix_to_list(IndexList);
    Ex=Expect;
  }
  return(Ex);
}

// setup for a distributed computation
def setup() {
#ifdef USE_MODULE
#else
    extern Ekn_plist;
    extern Ekn_IDL;
#endif

    if(getopt(report) != 1){
        if(type(getopt(nps)) >= 0){
    	    NPS = getopt(nps);
    	}else if(type(getopt(number_of_processes)) >= 0){
            NPS = getopt(number_of_processes);
        }else NPS = 1; //default

    	if(type(getopt(nprm)) >= 0){
            NPRM = getopt(nprm);
        }else if(type(getopt(number_of_primes)) >= 0){
            NPRM = getopt(number_of_primes);
        }else  NPRM = 10;  //default

    	if(type(getopt(fgp)) >= 0){
            FGP = getopt(fgp);
    	}else if(type(getopt(file_of_generated_primes)) >= 0){
            FGP = getopt(file_of_generated_primes);
    	}else  FGP = 0;//default

	if(type(getopt(minp)) >= 0){
	    MINP = getopt(minp);
	}else  MINP = 10^10;  //default

    	if(type(NPS) > 0){
            if(type(NPS) == 1){
	        if(type(Ekn_IDL) == 4)  shutdown(Ekn_IDL);
	    	Ekn_IDL = process(NPS|option_list=getopt());
	    }else  printf("Option nps is bad.\n");
    	}

    	if(type(NPRM) > 0){
            if(type(NPRM) == 7){
	        Ekn_plist = load_prime(NPRM);
	    }else if(type(NPRM) == 1){
	        if(type(MINP) == 1){
	            Ekn_plist = generate_p_set(MINP,NPRM);
		    if(type(FGP)==7){
		        save_prime(Ekn_plist,FGP);
		        printf("File name %a is written.\n",FGP);
		    }
 		    Ekn_plist = reverse(Ekn_plist);
	        }else    printf("Option minp is bad.\n");
	    }else{
	        printf("Option nprm is bad.\n");
	    }
        }
    }

    if(type(Ekn_IDL)==4)    printf("Number of processes = %a.\n",length(Ekn_IDL));
    else    printf("List of child processes is empty.\n");
    
    if(type(Ekn_plist)==4){
        printf("Number of primes = %a.\n",length(Ekn_plist));
	printf("Min of plist = %a.\n",Ekn_plist[0]);
    }else    printf("List of primes is empty\n");
    return;
}

def report(){
    setup(|report=1);
}

def cmle(U)
"It returns an approximate value of the conditional MLE for the table U. \
Options: see cBasistoEdiff.\
Example:   U=[[1,1,2,3],[1,3,1,1]]; M=gtt_ekn.cmle(U); M[2] is the probability. \
M[0] is U(input data), M[1] is maginal sums."
{
//    return(cmle_2(U|option_list=getopt()));
   return(cmle_3(U|option_list=getopt()));
}

def set_debug_level(Mode) {
  Ekn_debug=Mode;
}

def assert1(N) {
  if (N<2) return("N should be more than 1.");
  if (length(getopt()) == 0) printf("Try g_mat_fac_test_int: ");
  else printf("Try %a\n",getopt());
  T0=time();
  V1=gtt_ekn.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                       [[1,(1-1/N)/56],[1,1]]|option_list=getopt());
  V1_3=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                       [[1,(1-1/N)/56],[1,1]]|option_list=getopt());
  T1=time();
  printf("Timing =%a (CPU) + %a (GC) = %a (total), real time=%a\n",TT0=T1[0]-T0[0],TT1=T1[1]-T0[1],TT0+TT1,T1[3]-T0[3]);
  printf("\nTry g_mat_fac_test_plain: ");
  T0=time();
  V0=gtt_ekn.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                       [[1,(1-1/N)/56],[1,1]]|plain=1);
  V0_3=gtt_ekn3.gmvector([[36*N,13*N-1],[38*N-1,11*N]],
                       [[1,(1-1/N)/56],[1,1]]|plain=1);
  T1=time();
  printf("Timing =%a (CPU) + %a (GC) = %a (total)\n",TT0=T1[0]-T0[0],TT1=T1[1]-T0[1],TT0+TT1);
  printf("diff of both method and gtt_ekn and gtt_ekn3= \n");
  return([V0-V1,V0-V0_3,V1-V1_3]);
}

def assert2(N) {
  if (N<1) return("N should be more than 0.");
  if (N>1) printf("The second method will take time, e.g., N=2 ==> 1 min.\n");
  Marginal = [[130*N,170*N,353*N],[90*N,119*N,444*N]];
  P=[[17/100,1,10],[14/100,1,33/10],[1,1,1]];
  printf("Marginal=%a\nP=%a\n",Marginal,P);
  if (length(getopt()) == 0) printf("Try g_mat_fac_test_int: ");
  else printf("Try %a\n",getopt());
  T0=time();
  printf("Do V1_3=gtt_ekn3.expectation(Marginal=%a,P=%a|option_list=%a)\n",Marginal,P,getopt());
  V1_3=gtt_ekn3.expectation(Marginal,P|option_list=getopt())$
  T1=time();
  printf("gtt_ekn3: Timing (int) =%a (CPU) + %a (GC) = %a (total), real time=%a\n",TT0=T1[0]-T0[0],TT1=T1[1]-T0[1],TT0+TT1,T1[3]-T0[3]);
  T0=time();
  printf("Do V1=gtt_ekn.expectation(Marginal=%a,P=%a|option_list=%a)\n",Marginal,P,getopt());
  V1=gtt_ekn.expectation(Marginal,P|option_list=getopt())$
  T1=time();
  printf("gtt_ekn: Timing (int) =%a (CPU) + %a (GC) = %a (total), real time=%a\n",TT0=T1[0]-T0[0],TT1=T1[1]-T0[1],TT0+TT1,T1[3]-T0[3]);
  printf("\nTry g_mat_fac_test_plain: ");
  T0=time();
  printf("Do V0=gtt_ekn.expectation(Marginal=%a,P=%a|plain=1)\n",Marginal,P);
  V0=E=gtt_ekn.expectation(Marginal,P | plain=1);
  T1=time();
  printf("Timing (rational) =%a (CPU) + %a (GC) = %a (total)\n",TT0=T1[0]-T0[0],TT1=T1[1]-T0[1],TT0+TT1);
  T0=time();
  printf("Do V0_3=gtt_ekn3.expectation(Marginal=%a,P=%a|plain=1)\n",Marginal,P);
  V0_3=E=gtt_ekn3.expectation(Marginal,P | plain=1);
  T1=time();
  printf("Timing (rational) =%a (CPU) + %a (GC) = %a (total)\n",TT0=T1[0]-T0[0],TT1=T1[1]-T0[1],TT0+TT1);
  printf("diff of both method, gtt_ekn and gtt_ekn3 = \n");
  return([V0-V1,V1-V1_3,V0-V0_3]);
}

def assert3(R1,R2,Size) 
"assert3(R1,R2,Size | x=0, nps=3, crt=0); x=1 --> show windows of subprocesses. Example: gtt_ekn3.assert3(5,5,100 | x=1); inverval=m --> reduce every m steps of lin. trans."
{
  printf("Examples: assert3(R1=3,R2=5,Size=100|x=1); example3(5,5,100 | nps=32, crt=1);\n");
  printf("Note: load gtt_ekn3/ekn_eval-timing.rr to show timing data.\n");
  WorkDir="tmp-gtt_ekn3";
  printf("mkdir a temporary work directory %a\n",WorkDir);
  shell("mkdir -p "+WorkDir);
  P=prob1(R1,R2,Size);
  printf("Prob=%a\n",P);
  printf("Sequential method: ");
  G=gmvector(P[0],P[1] | option_list=getopt())$
  bsave(G,WorkDir+"/1.ab");
  printf("Done. Saved in 1.ab\n");

  printf("Parallel method: ");
  Nproc=3; if (type(getopt(nps))>0) Nproc=getopt(nps);
  if (type(getopt(x))>0) X=1; else X=0;
  printf("Number of process=%a, ",Nproc);
  gtt_ekn3.setup(|nps=Nproc,sub_progs=["gtt_ekn3.rr","gtt_ekn3/childprocess.rr"],x=X,
    fgp=WorkDir+"/p300.txt",nprm=300,minp=10^100);
  G2=gmvector(P[0],P[1] | option_list=getopt())$
  bsave(G,WorkDir+"/2.ab");
  printf("Done. Saved in 2.ab\n");
  printf("Diff (should be 0)=%a\n",base_flatten(number_eval(G-G2)));
}

/* bench mark test. */
def prob1(R1,R2,Size) {
  if (R1 > R2) return("Should be R1<=R2.");
  if (type(getopt(factor))>0) Factor=getopt(factor);
  else Factor=1;
  if (type(getopt(factor_row))>0) Factor_row=getopt(factor_row);
  else Factor_row=1;
  Col=newvect(R1);
  Row=newvect(R2);
  Sum2=0;
  for (I=0; I<R2; I++) {
    Row[I] = number_floor(Size*Factor_row*(I+1));
    Sum2 += Row[I];
  }
  Sum=0;
  for (I=0; I<R1-1; I++) {
    Col[I] = Size*Factor*(I+1);
    Sum += Col[I];
  }
  Col[R1-1] = Sum2-Sum;
  if (Row[R2-1] <= 0) printf("Error: negative marginal sum.\n");
  P=newmat(R1,R2);
  for (I=0; I<R1; I++) P[I][0] = 1;
  for (J=1; J<R2; J++) P[R1-1][J] =1;
  Q=2;
  for (I=0; I<R1-1; I++) {
    for (J=1; J<R2; J++) {
       P[I][J] = 1/(Q=pari(nextprime,Q)); Q++;
    }
  }
  return [[vtol(Col),vtol(Row)],matrix_matrix_to_list(P)];
}

def prob5x5(N) {
  Margin=[[4*N,4*N,4*N,4*N,4*N],[2*N,3*N,5*N,5*N,5*N]];
  P=[[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/37,1/41,1/43,1/47],[1,1,1,1,1]];
  return [Margin,P];
}

def set_ARule(Rule) { return(ARule=Rule); }
def set_XRule(Rule) { return(XRule=Rule); }

#ifdef USE_MODULE
endmodule$
gtt_ekn3.init_dm_bsplit()$
gtt_ekn3.init_bsplit(|minsize=16, levelmax=1)$
#else
#endif

import("gtt_ekn.rr")$  // for assert

end$