=================================================================== RCS file: /home/cvs/OpenXM/src/hgm/mh/src/jack-n.c,v retrieving revision 1.31 retrieving revision 1.46 diff -u -p -r1.31 -r1.46 --- OpenXM/src/hgm/mh/src/jack-n.c 2015/04/02 00:11:32 1.31 +++ OpenXM/src/hgm/mh/src/jack-n.c 2016/02/15 06:02:39 1.46 @@ -4,8 +4,10 @@ #include #include #include "sfile.h" + +#define VSTRING "%!version2.0" /* - $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.30 2015/03/24 07:49:06 takayama Exp $ + $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.45 2016/02/13 22:56:50 takayama Exp $ Ref: copied from this11/misc-2011/A1/wishart/Prog jack-n.c, translated from mh.rr or tk_jack.rr in the asir-contrib. License: LGPL Koev-Edelman for higher order derivatives. @@ -15,6 +17,11 @@ 2. Use the recurrence to obtain beta(). 3. Easier input data file format for mh-n.c Changelog: + 2016.02.15 log of Ef + 2016.02.09 unify 2F1 and 1F1. Parser. + 2016.02.04 Ef_type (exponential or scalar factor type) + 2016.02.01 ifdef C_2F1 ... + 2016.01.12 2F1 2014.03.15 http://fe.math.kobe-u.ac.jp/Movies/oxvh/2014-03-11-jack-n-c-automatic see also hgm/doc/ref.html, @s/2014/03/15-my-note-jack-automatic-order-F_A-casta.pdf. 2014.03.14, --automatic option. Output estimation data. 2012.02.21, porting to OpenXM/src/hgm @@ -55,8 +62,14 @@ static double Ef2; #define M_n0 3 /* used for tests. Must be equal to M_n */ #define M_m_MAX 200 #define M_nmx M_m_MAX /* maximal of M_n */ -#define A_LEN 1 /* (a_1) , (a_1, ..., a_p)*/ -#define B_LEN 1 /* (b_1) */ + +static int A_LEN=-1; /* (a_1, ..., a_p), A_LEN=p */ +static int B_LEN=-1; /* (b_1,..., b_q), B_LEN=q */ +static double *A_pFq=NULL; +static double *B_pFq=NULL; +static int Ef_type=1; /* 1F1 for Wishart */ +/* Ef_type=2; 2F1, Hashiguchi note (1) */ + static int Debug = 0; static int Alpha = 2; /* 2 implies the zonal polynomial */ static int *Darray = NULL; @@ -170,25 +183,18 @@ static int mysum(int L[]); static int plength(int P[]); static int plength_t(int P[]); static void ptrans(int P[M_nmx],int Pt[]); -static int test_ptrans(); static int huk(int K[],int I,int J); static int hdk(int K[],int I,int J); static double jjk(int K[]); static double ppoch(double A,int K[]); static double ppoch2(double A,double B,int K[]); static double mypower(double x,int n); -static double qk(int K[],double A[A_LEN],double B[B_LEN]); +static double qk(int K[],double A[],double B[]); static int bb(int N[],int K[],int M[],int I,int J); static double beta(int K[],int M[]); static int printp(int kappa[]); -static int printp2(int kappa[]); -static int test_beta(); static double q3_10(int K[],int M[],int SK); -static double q3_5(double A[],double B[],int K[],int I); -static void mtest4(); -static void mtest4b(); static int nk(int KK[]); -static int plength2(int P1[],int P2[]); static int myeq(int P1[],int P2[]); static int pListPartition(int M,int N); static int pListPartition2(int Less,int From,int To, int M); @@ -204,31 +210,94 @@ static int genDarray2(int M,int N); static int isHStrip(int Kap[],int Nu[]); static void hsExec_beta(void); static int genBeta(int Kap[]); -static int checkBeta1(); static int psublen(int Kap[],int Mu[]); static int genJack(int M,int N); -static int checkJack1(int M,int N); -static int checkJack2(int M,int N); -static int mtest1b(); static int imypower(int x,int n); static int usage(); static int setParamDefault(); static int next(struct SFILE *fp,char *s,char *msg); -static int gopen_file(void); static int setParam(char *fname); static int showParam(struct SFILE *fp,int fd); +#ifdef STANDALONE +static int showParam_v1(struct SFILE *fp,int fd); +#endif static double iv_factor(void); static double gammam(double a,int n); static double mypower(double a,int n); -double mh_t(double A[A_LEN],double B[B_LEN],int N,int M); +static double iv_factor_ef_type1(void); +static double iv_factor_ef_type2(void); +static double lgammam(double a,int n); +static double liv_factor_ef_type1(void); +static double liv_factor_ef_type2(void); + +static void setM_x(void); +static void setM_x_ef_type1(void); +static void setM_x_ef_type2(void); + +#ifdef STANDALONE +static int test_ptrans(); +static int printp2(int kappa[]); +static int test_beta(); +static void mtest4(); +static void mtest4b(); +static int plength2(int P1[],int P2[]); +static int checkBeta1(); +static int checkJack1(int M,int N); +static int checkJack2(int M,int N); +static int mtest1b(); +static double q3_5(double A[],double B[],int K[],int I); +#endif + +double mh_t(double A[],double B[],int N,int M); double mh_t2(int J); struct MH_RESULT *jk_main(int argc,char *argv[]); struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree); int jk_freeWorkArea(); int jk_initializeWorkArea(); +static void setA(double a[],int alen); /* set A_LEN and A_pFq */ +static void setB(double b[],int blen); +static void setA(double a[],int alen) { + int i; + if (alen < 0) { + if (A_pFq != NULL) myfree(A_pFq); + A_pFq=NULL; A_LEN=-1; + return; + } + if (alen == 0) { + A_LEN=0; return; + } + A_LEN=alen; + A_pFq = (double *)mymalloc(A_LEN*sizeof(double)); + if (a != NULL) { + for (i=0; i Q) { + for (I=0; I= M_kap[From], ..., M_kap[To], |(M_kap[From],...,M_kap[To])|<=M, */ static int pListPartition2(int Less,int From,int To, int M) { - int I,R; + int I; mh_check_intr(100); if (To < From) { (*M_pExec)(); return(0); } for (I=1; (I<=Less) && (I<=M) ; I++) { M_kap[From-1] = I; - R = pListPartition2(I,From+1,To,M-I); + pListPartition2(I,From+1,To,M-I); } return(1); } /* - partition に対してやる仕事をこちらへ書く. + Commands to do for each partition are given here. */ static void pExec_0() { if (Debug) { @@ -839,7 +926,7 @@ static int pListHS2(int From,int To,int Kap[]) { } static void hsExec_0() { - int i; + /* int i; */ if(Debug) {oxprintf("hsExec: "); printp(HS_mu); oxprintf("\n");} } @@ -1051,7 +1138,7 @@ static int genBeta(int Kap[]) { extern int M_beta_pt; extern int M_beta_kap[]; extern int P_pmn; - int I,J,N; + int I,N; if (Debug) {printp(Kap); oxprintf("<-Kappa, P_pmn=%d\n",P_pmn);} /* M_beta = newmat(2,P_pmn+1); */ M_beta_0 = (double *)mymalloc(sizeof(double)*(P_pmn+1)); @@ -1068,7 +1155,7 @@ static int genBeta(int Kap[]) { genBeta([2,2,0]); genBeta([2,1,1]); */ - +#ifdef STANDALONE static int checkBeta1() { int Kap[3] = {2,2,0}; int Kap2[3] = {2,1,0}; @@ -1099,7 +1186,7 @@ static int checkBeta1() { } return(0); } - +#endif /* def checkBeta2() { genDarray2(3,3); @@ -1153,7 +1240,7 @@ static int genJack(int M,int N) { int Pmn,I,J,K,L,Nv,H,P; int *Kap,*Mu; double Jack,Beta_km; - int Nk,JJ; + int Nk,JJ, two_to_I; if (Debug) oxprintf("genJack(%d,%d)\n",M,N); M_jack = (double **) mymalloc(sizeof(double *)*(N+2)); M_2n = imypower(2,N); @@ -1177,12 +1264,12 @@ static int genJack(int M,int N) { for (J=2; J= 2^I) aM_jack(I,J,K) = 0; + if (J >= two_to_I) aM_jack(I,J,K) = 0; /* J >= 2^I */ else aM_jack(I,J,K) = NAN; } } @@ -1254,7 +1341,7 @@ static int genJack(int M,int N) { } return(0); } - +#ifdef STANDALONE /* checkJack1(3,3) */ static int checkJack1(int M,int N) { @@ -1316,8 +1403,9 @@ static int checkJack2(int M,int N) { } /* main() { checkJack2(3,3); } */ +#endif -double mh_t(double A[A_LEN],double B[B_LEN],int N,int M) { +double mh_t(double A[],double B[],int N,int M) { double F,F2; extern int M_df; extern int P_pmn; @@ -1393,7 +1481,7 @@ double mh_t(double A[A_LEN],double B[B_LEN],int N,int oxprintf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g); oxprintfe("%%%%(stderr) serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m); oxprintfe("%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound); - oxprintfe("%%%%(stderr) F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g); + oxprintfe("%%%%(stderr) F=%lg,Ef=%lg,iv=|F*Ef|=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),iv,M_estimated_X0g,X0g); } M_mh_t_value=F; @@ -1414,6 +1502,7 @@ double mh_t2(int J) { return(F); } +#ifdef STANDALONE static int mtest1b() { double A[1] = {1.5}; double B[1] = {1.5+5}; @@ -1432,10 +1521,10 @@ static int mtest1b() { } /* main() { mtest1b(); }*/ +#endif - #define TEST 1 #ifndef TEST @@ -1456,7 +1545,7 @@ static int imypower(int x,int n) { } #ifdef STANDALONE2 -main(int argc,char *argv[]) { +int main(int argc,char *argv[]) { mh_exit(MH_RESET_EXIT); /* jk_main(argc,argv); oxprintf("second run.\n"); */ @@ -1489,11 +1578,11 @@ struct MH_RESULT *jk_main(int argc,char *argv[]) { } struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree) { - double *y0; - double x0,xn; - double ef; - int i,j,rank; - double a[1]; double b[1]; + // double *y0; + // double x0,xn; + // double ef; + + int i,j; // int i,j,rank; extern double M_x[]; extern double *Beta; extern int M_2n; @@ -1561,7 +1650,7 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a if (M_n > Mapprox) Mapprox=M_n; } /* Output by a file=stdout */ - ofp = mh_fopen("oxstdout","w",JK_byFile); + ofp = mh_fopen("stdout","w",JK_byFile); sprintf(swork,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp); sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp); @@ -1569,9 +1658,7 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a myerror("Mg must be equal to M_n\n"); mh_exit(-1); } if (Debug) showParam(NULL,1); - for (i=0; itest2.txt\n"); - oxprintfe(" ./hgm_w-n --idata test2.txt --gnuplotf test-g\n"); - oxprintfe(" gnuplot -persist s); + if (JK_byFile) oxprintfe("File %s is not found.\n",fname); mh_exit(-1); } - next(fp,s,"Mg(m)"); + /* set default initial values */ + Mg=-1; /* number of variables */ + Ng=(double *) mymalloc(sizeof(double)); *Ng=-1; /* *Ng is the degree of freedom 1F1 */ + X0g=0.1; /* evaluation point */ + Ef=1.0; /* exponential factor */ + Ef_type=1; + Hg=0.01; /* step size for RK */ + Dp = 1; /* sampling rate */ + Xng = 10.0; /* terminal point for RK */ + + /* Parser for the old style (version <2.0) input */ + version=next(fp,s,"Mg(m)"); + if (version == 2) goto myparse; sscanf(s,"%d",&Mg); rank = imypower(2,Mg); @@ -1727,7 +1804,6 @@ static int setParam(char *fname) { sscanf(s,"%lf",&(Beta[i])); } - Ng = (double *)mymalloc(sizeof(double)); next(fp,s,"%Ng= (freedom parameter n)"); sscanf(s,"%lf",Ng); @@ -1758,6 +1834,7 @@ static int setParam(char *fname) { sscanf(s,"%lf",&Xng); /* Reading the optional parameters */ + myparse: while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) { /* expect ID */ if (tk.type != MH_TOKEN_ID) { @@ -1823,24 +1900,212 @@ static int setParam(char *fname) { M_show_autosteps = tk.ival; continue; } + // Format: #p_pFq=2 1.5 3.2 + if (strcmp(s,"p_pFq")==0) { + if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); + } + if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); + } + setA(NULL,tk.ival); + for (i=0; i= 0)) { + if (A_LEN <1) setA(NULL,1); + if (B_LEN <1) setB(NULL,1); + A_pFq[0] = ((double)Mg+1.0)/2.0; + B_pFq[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */ + if (Debug) oxprintf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",A_pFq[0],B_pFq[0],M_n,Mapprox); + } mh_fclose(fp); return(0); } - -static int showParam(struct SFILE *fp,int fd) { +#ifdef STANDALONE +/* may remove */ +static int showParam_v1(struct SFILE *fp,int fd) { int rank,i; char swork[1024]; if (fd) { - fp = mh_fopen("oxstdout","w",1); + fp = mh_fopen("stdout","w",1); } rank = imypower(2,Mg); sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp); for (i=0; i= 0) { + sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp); + } sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp); for (i=0; i= 0) { + sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp); + } + sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp); + sprintf(swork,"%%Iv=\n"); mh_fputs(swork,fp); + for (i=0; i