=================================================================== RCS file: /home/cvs/OpenXM/src/hgm/mh/src/jack-n.c,v retrieving revision 1.24 retrieving revision 1.40 diff -u -p -r1.24 -r1.40 --- OpenXM/src/hgm/mh/src/jack-n.c 2014/03/16 00:00:46 1.24 +++ OpenXM/src/hgm/mh/src/jack-n.c 2016/02/04 06:56:05 1.40 @@ -5,7 +5,7 @@ #include #include "sfile.h" /* - $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.23 2014/03/15 11:23:58 takayama Exp $ + $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.39 2016/02/04 02:52:19 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 +15,9 @@ 2. Use the recurrence to obtain beta(). 3. Easier input data file format for mh-n.c Changelog: + 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 +58,26 @@ 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 */ +#ifdef C_2F1 +#define A_LEN 2 /* (a_1) , (a_1, ..., a_p)*/ +#define B_LEN 1 /* (b_1) */ +static int P_pFq=2; +static int Q_pFq=1; +#else #define A_LEN 1 /* (a_1) , (a_1, ..., a_p)*/ #define B_LEN 1 /* (b_1) */ +static int P_pFq=1; +static int Q_pFq=1; +#endif +static double A_pFq[A_LEN]; +static double B_pFq[B_LEN]; +#ifndef C_2F1 +static int Orig_1F1=1; +static int Ef_type=1; /* 1F1 for Wishart */ +#else +static int Orig_1F1=0; +static int Ef_type=2; /* 2F1, Hashiguchi note (1) */ +#endif static int Debug = 0; static int Alpha = 2; /* 2 implies the zonal polynomial */ static int *Darray = NULL; @@ -100,7 +121,7 @@ static double M_rel_error=0.0; /* relative errors */ If automatic == 1, then the series is reevaluated as long as t_success!=1 by increasing X0g (evaluation point) and M_m (approx degree); */ -static int M_automatic=0; +static int M_automatic=1; /* Estimated degree bound for series expansion. See mh_t */ static int M_m_estimated_approx_deg=0; /* Let F(i) be the approximation up to degree i. @@ -112,13 +133,13 @@ static double M_series_error; M_series_error < M_assigend_series_error (A) is required for the estimated_approx_deg. */ -static double M_assigned_series_error=0.00001; +static double M_assigned_series_error=M_ASSIGNED_SERIES_ERROR_DEFAULT; /* Let Ef be the exponential factor ( Ef=(4)/1F1 of [HNTT] ) If F(M_m)*Ef < x0value_min (B), the success=0 and X0g is increased. Note that minimal double is about 2e-308 */ -static double M_x0value_min=1e-30; +static double M_x0value_min=1e-60; /* estimated_X0g is the suggested value of X0g. */ @@ -138,6 +159,11 @@ static int M_mh_t_success=1; */ static double M_recommended_abserr; /* + recommended_relerr is the recommended value of the relative error + for the Runge-Kutta method.. + */ +static double M_recommended_relerr; +/* max of beta(i)*x/2 */ static double M_beta_i_x_o2_max; @@ -149,8 +175,11 @@ static double M_beta_i_beta_j_min; Value of matrix hg */ static double M_mh_t_value; +/* + Show the process of updating degree. + */ +int M_show_autosteps=1; - /* prototypes */ static void *mymalloc(int size); static int myfree(void *p); @@ -162,7 +191,6 @@ 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[]); @@ -173,14 +201,8 @@ static double qk(int K[],double A[A_LEN],double B[B_LE 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); @@ -196,24 +218,39 @@ 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); static double iv_factor(void); static double gammam(double a,int n); static double mypower(double a,int n); +static double iv_factor_ef_type1(void); +static double iv_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[A_LEN],double B[B_LEN],int N,int M); double mh_t2(int J); struct MH_RESULT *jk_main(int argc,char *argv[]); @@ -234,7 +271,7 @@ int jk_freeWorkArea() { if (M_beta_1) {myfree(M_beta_1); M_beta_1=NULL;} if (M_jack) { for (i=0; M_jack[i] != NULL; i++) { - if (Debug) printf("Free M_jack[%d]\n",i); + if (Debug) oxprintf("Free M_jack[%d]\n",i); myfree(M_jack[i]); M_jack[i] = NULL; } myfree(M_jack); M_jack=NULL; @@ -278,24 +315,24 @@ int jk_initializeWorkArea() { static void *mymalloc(int size) { void *p; - if (Debug) printf("mymalloc(%d)\n",size); + if (Debug) oxprintf("mymalloc(%d)\n",size); p = (void *)mh_malloc(size); if (p == NULL) { - fprintf(stderr,"No more memory.\n"); + oxprintfe("No more memory.\n"); mh_exit(-1); } return(p); } static int myfree(void *p) { - if (Debug) printf("myFree at %p\n",p); + if (Debug) oxprintf("myFree at %p\n",p); return(mh_free(p)); } -static int myerror(char *s) { fprintf(stderr,"%s: type in control-C\n",s); getchar(); getchar(); return(0);} +static int myerror(char *s) { oxprintfe("%s: type in control-C\n",s); getchar(); getchar(); return(0);} static double jack1(int K) { double F; extern int Alpha; - int I,J,L,II,JJ,N; + int J,II,JJ,N; /* int I,J,L,II,JJ,N; */ N = 1; if (K == 0) return((double)1); F = xval(1,K); @@ -308,7 +345,7 @@ static double jack1(int K) { static double jack1diff(int K) { double F; extern int Alpha; - int I,J,S,L,II,JJ,N; + int J,II,JJ,N; /* int I,J,S,L,II,JJ,N; */ N = 1; if (K == 0) return((double)1); F = K*xval(1,K-1); @@ -321,7 +358,6 @@ static double jack1diff(int K) { static double xval(int ii,int p) { /* x_i^p */ extern double M_x[]; - double F; int i,j; static int init=0; if (JK_deallocate) { init=0; return(0.0);} @@ -341,7 +377,7 @@ static double xval(int ii,int p) { /* x_i^p */ if (p > M_m_MAX-2) myerror("xval, p is too large."); if (p < 0) { myerror("xval, p is negative."); - printf("ii=%d, p=%d\n",ii,p); + oxprintf("ii=%d, p=%d\n",ii,p); mh_exit(-1); } return(Xarray[ii-1][p]); @@ -390,6 +426,7 @@ static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] } } +#ifdef STANDALONE static int test_ptrans() { extern int M_m; int p[M_n0]={5,3,2}; @@ -397,11 +434,11 @@ static int test_ptrans() { int i; M_m = 10; ptrans(p,pt); - if (Debug) {for (i=0; i<10; i++) printf("%d,",pt[i]); printf("\n");} + if (Debug) {for (i=0; i<10; i++) oxprintf("%d,",pt[i]); oxprintf("\n");} return(0); } +#endif - /* upper hook length h_kappa^*(K) @@ -411,7 +448,7 @@ static int huk(int K[],int I,int J) { int Kp[M_m_MAX]; int A,H; A=Alpha; - /*printf("h^k(%a,%a,%a,%a)\n",K,I,J,A);*/ + /*oxprintf("h^k(%a,%a,%a,%a)\n",K,I,J,A);*/ ptrans(K,Kp); H=Kp[J-1]-I+A*(K[I-1]-J+1); return(H); @@ -426,7 +463,7 @@ static int hdk(int K[],int I,int J) { int Kp[M_m_MAX]; int A,H; A = Alpha; - /*printf("h_k(%a,%a,%a,%a)\n",K,I,J,A);*/ + /*oxprintf("h_k(%a,%a,%a,%a)\n",K,I,J,A);*/ ptrans(K,Kp); H=Kp[J-1]-I+1+A*(K[I-1]-J); return(H); @@ -436,9 +473,9 @@ static int hdk(int K[],int I,int J) { */ static double jjk(int K[]) { extern int Alpha; - int A,L,I,J; + int L,I,J; double V; - A=Alpha; + V=1; L=plength(K); 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) { - printf("M_kap="); + oxprintf("M_kap="); printp(M_kap); - printf("\n"); + oxprintf("\n"); } } @@ -831,8 +884,8 @@ static int pListHS2(int From,int To,int Kap[]) { } static void hsExec_0() { - int i; - if(Debug) {printf("hsExec: "); printp(HS_mu); printf("\n");} + /* int i; */ + if(Debug) {oxprintf("hsExec: "); printp(HS_mu); oxprintf("\n");} } /* @@ -870,10 +923,10 @@ static int pmn(int M,int N) { } P_pmn=S; if (Debug) { - printf("P_pmn=%d\n",P_pmn); + oxprintf("P_pmn=%d\n",P_pmn); for (i=0; i<=Min_m_n; i++) { - for (j=0; j<=M; j++) printf("%d,",aP_pki(i,j)); - printf("\n"); + for (j=0; j<=M; j++) oxprintf("%d,",aP_pki(i,j)); + oxprintf("\n"); } } myfree(P_pki); P_pki=NULL; @@ -881,7 +934,7 @@ static int pmn(int M,int N) { } /* - main() {pmn(4,3); printf("P_pmn=%d\n",P_pmn);} + main() {pmn(4,3); oxprintf("P_pmn=%d\n",P_pmn);} */ static int *cloneP(int a[]) { @@ -922,7 +975,7 @@ static int genDarray2(int M,int N) { M_m = M; Pmn = pmn(M,N)+1; - if (Debug) printf("Degree M = %d, N of vars N = %d, Pmn+1=%d\n",M,N,Pmn); + if (Debug) oxprintf("Degree M = %d, N of vars N = %d, Pmn+1=%d\n",M,N,Pmn); Darray=(int *) mymalloc(sizeof(int)*Pmn); for (i=0; i= M_n) { - if (Debug) {fprintf(stderr,"Ksize >= M_n\n");} + if (Debug) {oxprintfe("Ksize >= M_n\n");} continue; } for (i=0; i= 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; } } @@ -1191,7 +1244,7 @@ static int genJack(int M,int N) { for (J=1; J0) partial_sum[i] += partial_sum[i-1]; - serror = myabs((partial_sum[i]-partial_sum[i-1])/partial_sum[i-1]); + if (i>0) serror = myabs((partial_sum[i]-partial_sum[i-1])/partial_sum[i-1]); if ((i>0)&&(M_m_estimated_approx_deg < 0)&&(serror= argc) { fprintf(stderr,"Option argument is not given.\n"); return(NULL); }} +#define inci(i) { i++; if (i >= argc) { oxprintfe("Option argument is not given.\n"); return(NULL); }} static int imypower(int x,int n) { int i; @@ -1444,10 +1503,10 @@ 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); - printf("second run.\n"); */ + oxprintf("second run.\n"); */ jk_main(argc,argv); return(0); } @@ -1477,11 +1536,12 @@ 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; + double a[A_LEN]; double b[B_LEN]; extern double M_x[]; extern double *Beta; extern int M_2n; @@ -1511,7 +1571,7 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a }else if (strcmp(argv[i],"--help")==0) { usage(); return(0); }else if (strcmp(argv[i],"--bystring")==0) { - if (idata) {fprintf(stderr,"--bystring must come before --idata option.\n"); mh_exit(-1);} + if (idata) {oxprintfe("--bystring must come before --idata option.\n"); mh_exit(-1);} JK_byFile = 0; }else if (strcmp(argv[i],"--automatic")==0) { inci(i); /* ignore, in this function */ @@ -1522,7 +1582,7 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a inci(i); sscanf(argv[i],"%lg",&M_x0value_min); }else { - fprintf(stderr,"Unknown option %s\n",argv[i]); + oxprintfe("Unknown option %s\n",argv[i]); usage(); return(NULL); } @@ -1544,6 +1604,10 @@ struct MH_RESULT *jk_main2(int argc,char *argv[],int a ans->recommended_abserr = 1.0e-10; } else ans = NULL; + if (M_automatic) { + /* Differentiation can be M_m in the bit pattern in the M_n variable case.*/ + if (M_n > Mapprox) Mapprox=M_n; + } /* Output by a file=stdout */ ofp = mh_fopen("stdout","w",JK_byFile); @@ -1553,9 +1617,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; iseries_error = M_series_error; ans->recommended_abserr = M_recommended_abserr; } - if (Debug) printf("jk_freeWorkArea() starts\n"); + if (Debug) oxprintf("jk_freeWorkArea() starts\n"); jk_freeWorkArea(); - if (Debug) printf("jk_freeWorkArea() has finished.\n"); + if (Debug) oxprintf("jk_freeWorkArea() has finished.\n"); return(ans); } static int usage() { - fprintf(stderr,"Usages:\n"); - fprintf(stderr,"hgm_jack-n [--idata input_data_file --x0 x0 --degree approxm]\n"); - fprintf(stderr," [--automatic n --assigned_series_error e --x0value_min e2]\n"); - fprintf(stderr,"\nThe command hgm_jack-n [options] generates an input for hgm_w-n, Pr({y | ytest2.txt\n"); - fprintf(stderr," ./hgm_w-n --idata test2.txt --gnuplotf test-g\n"); - fprintf(stderr," gnuplot -persist test2.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",fp->s); mh_exit(-1); } next(fp,s,"Mg(m)"); @@ -1747,59 +1834,119 @@ static int setParam(char *fname) { while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) { /* expect ID */ if (tk.type != MH_TOKEN_ID) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } if (strcmp(s,"automatic")==0) { if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } M_automatic = tk.ival; continue; } if (strcmp(s,"assigned_series_error")==0) { if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } M_assigned_series_error = tk.dval; continue; } if (strcmp(s,"x0value_min")==0) { if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } M_x0value_min = tk.dval; continue; } if ((strcmp(s,"Mapprox")==0) || (strcmp(s,"degree")==0)) { if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } Mapprox = tk.ival; continue; } if (strcmp(s,"X0g_bound")==0) { if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) { - fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1); + oxprintfe("Syntax error at %s\n",s); mh_exit(-1); } M_X0g_bound = tk.dval; continue; } - fprintf(stderr,"Unknown ID at %s\n",s); mh_exit(-1); + if (strcmp(s,"show_autosteps")==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); + } + M_show_autosteps = tk.ival; + continue; + } + // Format: #p_pFq=2 1.5 3.2 + if (strcmp(s,"p_pFq")==0) { + Orig_1F1=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); + } + P_pFq = tk.ival; + for (i=0; i