Annotation of OpenXM/src/hgm/mh/src/jack-n.c, Revision 1.28
1.1 takayama 1: #include <stdio.h>
2: #include <stdlib.h>
3: #define _ISOC99_SOURCE
4: #include <math.h>
5: #include <string.h>
6: #include "sfile.h"
7: /*
1.28 ! takayama 8: $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.27 2014/03/20 09:37:16 takayama Exp $
1.12 takayama 9: Ref: copied from this11/misc-2011/A1/wishart/Prog
10: jack-n.c, translated from mh.rr or tk_jack.rr in the asir-contrib. License: LGPL
11: Koev-Edelman for higher order derivatives.
12: cf. 20-my-note-mh-jack-n.pdf, /Movies/oxvh/priv-j1/2011-12-14-ke-mh-jack.mov
13: Todo:
14: 1. Cash the transposes of partitions.
15: 2. Use the recurrence to obtain beta().
16: 3. Easier input data file format for mh-n.c
17: Changelog:
1.22 takayama 18: 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.
1.18 takayama 19: 2014.03.14, --automatic option. Output estimation data.
1.12 takayama 20: 2012.02.21, porting to OpenXM/src/hgm
21: 2011.12.22, --table option, which is experimental.
22: 2011.12.19, bug fix. jjk() should return double. It can become more than max int.
23: 2011.12.15, mh.r --> jack-n.c
1.1 takayama 24: */
25:
1.3 takayama 26: /****** from mh-n.c *****/
27: static int JK_byFile=1;
1.6 takayama 28: static int JK_deallocate=0;
1.3 takayama 29: #define M_n_default 3
1.4 takayama 30: #define Sample_default 1
1.3 takayama 31: static int M_n=0;
1.1 takayama 32: /* global variables. They are set in setParam() */
1.2 takayama 33: static int Mg; /* n */
34: static int Mapprox; /* m, approximation degree */
35: static double *Beta; /* beta[0], ..., beta[m-1] */
36: static double *Ng; /* freedom n. c=(m+1)/2+n/2; Note that it is a pointer */
37: static double X0g; /* initial point */
38: static double *Iv; /* Initial values of mhg sorted by mhbase() in rd.rr at beta*x0 */
39: static double Ef; /* exponential factor at beta*x0 */
40: static double Hg; /* step size of rk defined in rk.c */
41: static int Dp; /* Data sampling period */
42: static double Xng=0.0; /* the last point */
1.4 takayama 43: static int Sample = Sample_default;
1.1 takayama 44:
45: /* for sample inputs */
1.2 takayama 46: static double *Iv2;
47: static double Ef2;
1.1 takayama 48:
49: #ifdef NAN
50: #else
1.2 takayama 51: #define NAN 3.40282e+38 /* for old 32 bit machines. Todo, configure */
1.1 takayama 52: #endif
53:
54: /* #define M_n 3 defined in the Makefile */ /* number of variables */
55: #define M_n0 3 /* used for tests. Must be equal to M_n */
56: #define M_m_MAX 200
1.3 takayama 57: #define M_nmx M_m_MAX /* maximal of M_n */
1.1 takayama 58: #define A_LEN 1 /* (a_1) , (a_1, ..., a_p)*/
59: #define B_LEN 1 /* (b_1) */
1.9 takayama 60: static int Debug = 0;
1.5 takayama 61: static int Alpha = 2; /* 2 implies the zonal polynomial */
1.2 takayama 62: static int *Darray = NULL;
63: static int **Parray = NULL; /* array of partitions of size M_n */
64: static int *ParraySize = NULL; /* length of each partitions */
1.3 takayama 65: static int M_kap[M_nmx];
1.2 takayama 66: static int M_m=M_m_MAX-2; /* | | <= M_m, bug do check of M_m <=M_m_MAX-2 */
1.1 takayama 67: void (*M_pExec)(void);
1.3 takayama 68: static int HS_mu[M_nmx];
69: static int HS_n=M_nmx; /* It is initialized to M_n in jk_main */
1.1 takayama 70: void (*HS_hsExec)(void);
1.3 takayama 71: static double M_x[M_nmx];
1.1 takayama 72:
73: /* They are used in pmn */
1.2 takayama 74: static int *P_pki=NULL;
75: static int P_pmn=0;
1.1 takayama 76:
77: /* It is used genDarray2(), list partitions... */
1.2 takayama 78: static int DR_parray=0;
1.1 takayama 79:
80: /* Used in genBeta() and enumeration of horizontal strip. */
1.2 takayama 81: static double *M_beta_0=NULL; /* M_beta[0][*] value of beta_{kappa,mu}, M_beta[1][*] N_mu */
82: static int *M_beta_1=NULL;
83: static int M_beta_pt=0;
1.3 takayama 84: static int M_beta_kap[M_nmx];
1.2 takayama 85: static int UseTable = 0;
86:
87: static double **M_jack;
88: static int M_df=1; /* Compute differentials? */
89: static int M_2n=0; /* 2^N */
1.1 takayama 90:
1.3 takayama 91: static double Xarray[M_nmx][M_m_MAX];
1.1 takayama 92: /* (x_1, ..., x_n) */
93: /* Xarray[i][0] x_{i+1}^0, Xarray[i][1], x_{i+1}^1, ... */
94:
1.2 takayama 95: static double *M_qk=NULL; /* saves pochhammerb */
96: static double M_rel_error=0.0; /* relative errors */
1.1 takayama 97:
1.16 takayama 98: /* For automatic degree and X0g setting. */
1.19 takayama 99: /*
100: If automatic == 1, then the series is reevaluated as long as t_success!=1
101: by increasing X0g (evaluation point) and M_m (approx degree);
102: */
1.26 takayama 103: static int M_automatic=1;
1.19 takayama 104: /* Estimated degree bound for series expansion. See mh_t */
1.16 takayama 105: static int M_m_estimated_approx_deg=0;
1.19 takayama 106: /* Let F(i) be the approximation up to degree i.
107: The i-th series_error is defined
108: by |(F(i)-F(i-1))/F(i-1)|.
109: */
110: static double M_series_error;
111: /*
112: M_series_error < M_assigend_series_error (A) is required for the
113: estimated_approx_deg.
114: */
1.26 takayama 115: static double M_assigned_series_error=M_ASSIGNED_SERIES_ERROR_DEFAULT;
1.19 takayama 116: /*
117: Let Ef be the exponential factor ( Ef=(4)/1F1 of [HNTT] )
118: If F(M_m)*Ef < x0value_min (B), the success=0 and X0g is increased.
119: Note that minimal double is about 2e-308
120: */
1.26 takayama 121: static double M_x0value_min=1e-60;
1.19 takayama 122: /*
123: estimated_X0g is the suggested value of X0g.
124: */
1.16 takayama 125: static double M_estimated_X0g=0.0;
1.19 takayama 126: /*
127: X0g should be less than M_X0g_bound.
128: */
129: static double M_X0g_bound = 1e+100;
130: /*
131: success is set to 1 when (A) and (B) are satisfied.
132: */
1.16 takayama 133: static int M_mh_t_success=1;
1.19 takayama 134: /*
135: recommended_abserr is the recommended value of the absolute error
136: for the Runge-Kutta method. It is defined as
137: assigend_series_error(standing for significant digits)*Ig[0](initial value)
138: */
139: static double M_recommended_abserr;
140: /*
1.27 takayama 141: recommended_relerr is the recommended value of the relative error
142: for the Runge-Kutta method..
143: */
144: static double M_recommended_relerr;
145: /*
1.19 takayama 146: max of beta(i)*x/2
147: */
148: static double M_beta_i_x_o2_max;
149: /*
150: minimum of |beta_i-beta_j|
151: */
152: static double M_beta_i_beta_j_min;
153: /*
154: Value of matrix hg
155: */
156: static double M_mh_t_value;
1.25 takayama 157: /*
158: Show the process of updating degree.
159: */
160: int M_show_autosteps=1;
1.16 takayama 161:
1.2 takayama 162: /* prototypes */
163: static void *mymalloc(int size);
1.10 takayama 164: static int myfree(void *p);
165: static int myerror(char *s);
1.2 takayama 166: static double jack1(int K);
167: static double jack1diff(int k);
168: static double xval(int i,int p); /* x_i^p */
169: static int mysum(int L[]);
170: static int plength(int P[]);
171: static int plength_t(int P[]);
1.3 takayama 172: static void ptrans(int P[M_nmx],int Pt[]);
1.10 takayama 173: static int test_ptrans();
1.2 takayama 174: static int huk(int K[],int I,int J);
175: static int hdk(int K[],int I,int J);
176: static double jjk(int K[]);
177: static double ppoch(double A,int K[]);
178: static double ppoch2(double A,double B,int K[]);
179: static double mypower(double x,int n);
180: static double qk(int K[],double A[A_LEN],double B[B_LEN]);
181: static int bb(int N[],int K[],int M[],int I,int J);
182: static double beta(int K[],int M[]);
1.10 takayama 183: static int printp(int kappa[]);
184: static int printp2(int kappa[]);
185: static int test_beta();
1.2 takayama 186: static double q3_10(int K[],int M[],int SK);
187: static double q3_5(double A[],double B[],int K[],int I);
188: static void mtest4();
189: static void mtest4b();
190: static int nk(int KK[]);
191: static int plength2(int P1[],int P2[]);
192: static int myeq(int P1[],int P2[]);
193: static int pListPartition(int M,int N);
194: static int pListPartition2(int Less,int From,int To, int M);
195: static void pExec_0();
196: static int pListHS(int Kap[],int N);
197: static int pListHS2(int From,int To,int Kap[]);
198: static void hsExec_0();
199: static int pmn(int M,int N);
200: static int *cloneP(int a[]);
1.10 takayama 201: static int copyP(int p[],int a[]);
1.2 takayama 202: static void pExec_darray(void);
1.10 takayama 203: static int genDarray2(int M,int N);
204: static int isHStrip(int Kap[],int Nu[]);
1.2 takayama 205: static void hsExec_beta(void);
1.10 takayama 206: static int genBeta(int Kap[]);
207: static int checkBeta1();
1.2 takayama 208: static int psublen(int Kap[],int Mu[]);
1.10 takayama 209: static int genJack(int M,int N);
210: static int checkJack1(int M,int N);
211: static int checkJack2(int M,int N);
212: static int mtest1b();
1.2 takayama 213:
214: static int imypower(int x,int n);
1.10 takayama 215: static int usage();
216: static int setParamDefault();
217: static int next(struct SFILE *fp,char *s,char *msg);
1.2 takayama 218: static int gopen_file(void);
219: static int setParam(char *fname);
1.8 takayama 220: static int showParam(struct SFILE *fp,int fd);
1.2 takayama 221: static double iv_factor(void);
222: static double gammam(double a,int n);
223: static double mypower(double a,int n);
224:
225: double mh_t(double A[A_LEN],double B[B_LEN],int N,int M);
226: double mh_t2(int J);
1.3 takayama 227: struct MH_RESULT *jk_main(int argc,char *argv[]);
1.16 takayama 228: struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree);
1.3 takayama 229: int jk_freeWorkArea();
230: int jk_initializeWorkArea();
231:
232: int jk_freeWorkArea() {
233: /* bug, p in the cloneP will not be deallocated.
1.12 takayama 234: Nk in genDarray2 will not be deallocated.
235: */
1.3 takayama 236: int i;
1.6 takayama 237: JK_deallocate=1;
1.3 takayama 238: if (Darray) {myfree(Darray); Darray=NULL;}
239: if (Parray) {myfree(Parray); Parray=NULL;}
240: if (ParraySize) {myfree(ParraySize); ParraySize=NULL;}
241: if (M_beta_0) {myfree(M_beta_0); M_beta_0=NULL;}
242: if (M_beta_1) {myfree(M_beta_1); M_beta_1=NULL;}
243: if (M_jack) {
1.12 takayama 244: for (i=0; M_jack[i] != NULL; i++) {
245: if (Debug) printf("Free M_jack[%d]\n",i);
246: myfree(M_jack[i]); M_jack[i] = NULL;
247: }
248: myfree(M_jack); M_jack=NULL;
1.3 takayama 249: }
250: if (M_qk) {myfree(M_qk); M_qk=NULL;}
251: if (P_pki) {myfree(P_pki); P_pki=NULL;}
1.6 takayama 252: JK_deallocate=0;
1.13 takayama 253: return(0);
1.3 takayama 254: }
255: int jk_initializeWorkArea() {
256: int i,j;
1.6 takayama 257: JK_deallocate=1;
258: xval(0,0);
259: JK_deallocate=0;
1.3 takayama 260: Darray=NULL;
261: Parray=NULL;
262: ParraySize=NULL;
263: M_beta_0=NULL;
264: M_beta_1=NULL;
265: M_jack=NULL;
266: M_qk=NULL;
267: for (i=0; i<M_nmx; i++) M_kap[i]=HS_mu[i]=0;
268: for (i=0; i<M_nmx; i++) M_x[i]=0;
269: for (i=0; i<M_nmx; i++) for (j=0; j<M_m_MAX; j++) Xarray[i][j]=0;
1.5 takayama 270: for (i=0; i<M_nmx; i++) M_beta_kap[i]=0;
1.3 takayama 271: M_m=M_m_MAX-2;
272: Alpha = 2;
273: HS_n=M_nmx;
274: P_pki=NULL;
275: P_pmn=0;
276: DR_parray=0;
277: M_beta_pt=0;
278: M_df=1;
279: M_2n=0;
280: M_rel_error=0.0;
1.4 takayama 281: Sample = Sample_default;
1.5 takayama 282: Xng=0.0;
283: M_n=0;
1.13 takayama 284: return(0);
1.3 takayama 285: }
1.2 takayama 286:
287: static void *mymalloc(int size) {
1.1 takayama 288: void *p;
289: if (Debug) printf("mymalloc(%d)\n",size);
1.8 takayama 290: p = (void *)mh_malloc(size);
1.1 takayama 291: if (p == NULL) {
292: fprintf(stderr,"No more memory.\n");
1.3 takayama 293: mh_exit(-1);
1.1 takayama 294: }
295: return(p);
296: }
1.10 takayama 297: static int myfree(void *p) {
1.8 takayama 298: if (Debug) printf("myFree at %p\n",p);
1.13 takayama 299: return(mh_free(p));
1.7 takayama 300: }
1.13 takayama 301: static int myerror(char *s) { fprintf(stderr,"%s: type in control-C\n",s); getchar(); getchar(); return(0);}
1.1 takayama 302:
1.2 takayama 303: static double jack1(int K) {
1.1 takayama 304: double F;
305: extern int Alpha;
306: int I,J,L,II,JJ,N;
307: N = 1;
308: if (K == 0) return((double)1);
309: F = xval(1,K);
310: for (J=0; J<K; J++) {
311: II = 1; JJ = J+1;
312: F *= (N-(II-1)+Alpha*(JJ-1));
313: }
314: return(F);
315: }
1.2 takayama 316: static double jack1diff(int K) {
1.1 takayama 317: double F;
318: extern int Alpha;
319: int I,J,S,L,II,JJ,N;
320: N = 1;
321: if (K == 0) return((double)1);
322: F = K*xval(1,K-1);
323: for (J=0; J<K; J++) {
324: II = 1; JJ = J+1;
325: F *= (N-(II-1)+Alpha*(JJ-1));
326: }
327: return(F);
328: }
329:
1.2 takayama 330: static double xval(int ii,int p) { /* x_i^p */
1.1 takayama 331: extern double M_x[];
332: double F;
333: int i,j;
1.10 takayama 334: static int init=0;
1.6 takayama 335: if (JK_deallocate) { init=0; return(0.0);}
1.1 takayama 336: if (!init) {
337: for (i=1; i<=M_n; i++) {
338: for (j=0; j<M_m_MAX; j++) {
1.12 takayama 339: if (j != 0) {
340: Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1];
341: }else{
342: Xarray[i-1][0] = 1;
343: }
1.1 takayama 344: }
345: }
346: init = 1;
347: }
348: if (ii < 1) myerror("xval, index out of bound.");
349: if (p > M_m_MAX-2) myerror("xval, p is too large.");
350: if (p < 0) {
351: myerror("xval, p is negative.");
352: printf("ii=%d, p=%d\n",ii,p);
1.3 takayama 353: mh_exit(-1);
1.1 takayama 354: }
355: return(Xarray[ii-1][p]);
356: }
357:
1.2 takayama 358: static int mysum(int L[]) {
1.1 takayama 359: int S,I,N;
360: N=M_n;
361: S=0;
362: for (I=0; I<N; I++) S += L[I];
363: return(S);
364: }
365:
366: /*
1.12 takayama 367: (3,2,2,0,0) --> 3
1.1 takayama 368: */
1.2 takayama 369: static int plength(int P[]) {
1.1 takayama 370: int I;
371: for (I=0; I<M_n; I++) {
372: if (P[I] == 0) return(I);
373: }
374: return(M_n);
375: }
376: /* plength for transpose */
1.2 takayama 377: static int plength_t(int P[]) {
1.1 takayama 378: int I;
379: for (I=0; I<M_m; I++) {
380: if (P[I] == 0) return(I);
381: }
382: return(M_m);
383: }
384:
385: /*
1.12 takayama 386: ptrans(P) returns Pt
1.1 takayama 387: */
1.3 takayama 388: static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] */
1.1 takayama 389: extern int M_m;
390: int i,j,len;
1.3 takayama 391: int p[M_nmx];
1.1 takayama 392: for (i=0; i<M_n; i++) p[i] = P[i];
393: for (i=0; i<M_m+1; i++) Pt[i] = 0;
394: for (i=0; i<M_m; i++) {
395: len=plength(p); Pt[i] = len;
396: if (len == 0) return;
397: for (j=0; j<len; j++) p[j] -= 1;
398: }
399: }
400:
1.10 takayama 401: static int test_ptrans() {
1.1 takayama 402: extern int M_m;
403: int p[M_n0]={5,3,2};
404: int pt[10];
405: int i;
406: M_m = 10;
407: ptrans(p,pt);
408: if (Debug) {for (i=0; i<10; i++) printf("%d,",pt[i]); printf("\n");}
1.13 takayama 409: return(0);
1.1 takayama 410: }
411:
412:
413: /*
414: upper hook length
415: h_kappa^*(K)
416: */
1.2 takayama 417: static int huk(int K[],int I,int J) {
1.1 takayama 418: extern int Alpha;
419: int Kp[M_m_MAX];
420: int A,H;
421: A=Alpha;
422: /*printf("h^k(%a,%a,%a,%a)\n",K,I,J,A);*/
423: ptrans(K,Kp);
424: H=Kp[J-1]-I+A*(K[I-1]-J+1);
425: return(H);
426: }
427:
428: /*
429: lower hook length
430: h^kappa_*(K)
431: */
1.2 takayama 432: static int hdk(int K[],int I,int J) {
1.1 takayama 433: extern int Alpha;
434: int Kp[M_m_MAX];
435: int A,H;
436: A = Alpha;
437: /*printf("h_k(%a,%a,%a,%a)\n",K,I,J,A);*/
438: ptrans(K,Kp);
439: H=Kp[J-1]-I+1+A*(K[I-1]-J);
440: return(H);
441: }
442: /*
443: j_kappa. cf. Stanley.
444: */
1.2 takayama 445: static double jjk(int K[]) {
1.1 takayama 446: extern int Alpha;
447: int A,L,I,J;
448: double V;
449: A=Alpha;
450: V=1;
451: L=plength(K);
452: for (I=0; I<L; I++) {
453: for (J=0; J<K[I]; J++) {
454: V *= huk(K,I+1,J+1)*hdk(K,I+1,J+1);
455: }
456: }
457: if (Debug) {printp(K); printf("<--K, jjk=%lg\n",V);}
458: return(V);
459: }
460: /*
461: (a)_kappa^\alpha, Pochhammer symbol
462: Note that ppoch(a,[n]) = (a)_n, Alpha=2
463: */
1.2 takayama 464: static double ppoch(double A,int K[]) {
1.1 takayama 465: extern int Alpha;
466: double V;
467: int L,I,J,II,JJ;
468: V = 1;
469: L=plength(K);
470: for (I=0; I<L; I++) {
471: for (J=0; J<K[I]; J++) {
472: II = I+1; JJ = J+1;
473: V *= (A-((double)(II-1))/((double)Alpha)+JJ-1);
474: }
475: }
476: return(V);
477: }
1.2 takayama 478: static double ppoch2(double A,double B,int K[]) {
1.1 takayama 479: extern int Alpha;
480: double V;
481: int L,I,J,II,JJ;
482: V = 1;
483: L=plength(K);
484: for (I=0; I<L; I++) {
485: for (J=0; J<K[I]; J++) {
486: II = I+1; JJ = J+1;
487: V *= (A-((double)(II-1))/((double)Alpha)+JJ-1);
488: V /= (B-((double)(II-1))/((double)Alpha)+JJ-1);
489: }
490: }
491: return(V);
492: }
1.2 takayama 493: static double mypower(double x,int n) {
1.1 takayama 494: int i;
495: double v;
496: if (n < 0) return(1/mypower(x,-n));
497: v = 1;
498: for (i=0; i<n; i++) v *= x;
499: return(v);
500: }
501: /* Q_kappa
1.12 takayama 502: */
1.2 takayama 503: static double qk(int K[],double A[A_LEN],double B[B_LEN]) {
1.1 takayama 504: extern int Alpha;
505: int P,Q,I;
506: double V;
507: P = A_LEN;
508: Q = B_LEN;
509: V = mypower((double) Alpha,mysum(K))/jjk(K);
1.2 takayama 510: /* to reduce numerical errors, temporary. */
1.1 takayama 511: if (P == Q) {
512: for (I=0; I<P; I++) V = V*ppoch2(A[I],B[I],K);
513: }
514: return(V);
515:
516: for (I=0; I<P; I++) V = V*ppoch(A[I],K);
517: for (I=0; I<Q; I++) V = V/ppoch(B[I],K);
518: return(V);
519: }
520:
521: /*
1.12 takayama 522: B^nu_{kappa,mu}(i,j)
523: bb(N,K,M,I,J)
1.1 takayama 524: */
1.2 takayama 525: static int bb(int N[],int K[],int M[],int I,int J) {
1.1 takayama 526: int Kp[M_m_MAX]; int Mp[M_m_MAX];
527: ptrans(K,Kp);
528: ptrans(M,Mp);
529:
530: /*
1.12 takayama 531: printp(K); printf("K<--, "); printp2(Kp); printf("<--Kp\n");
532: printp(M); printf("M<--, "); printp2(Mp); printf("<--Mp\n");
1.1 takayama 533: */
534:
535: if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J));
536: if (Kp[J-1] == Mp[J-1]) return(huk(N,I,J));
537: else return(hdk(N,I,J));
538: }
539: /*
540: beta_{kappa,mu}
541: beta(K,M)
542: */
1.2 takayama 543: static double beta(int K[],int M[]) {
1.1 takayama 544: double V;
545: int L,I,J,II,JJ;
546: V = 1;
547:
548: L=plength(K);
549: for (I=0; I<L; I++) {
550: for (J=0; J<K[I]; J++) {
551: II = I+1; JJ = J+1;
552: V *= (double)bb(K,K,M,II,JJ);
553: /* printf("[%d,%d,%lf]\n",I,J,V); */
554: }
555: }
556:
557: L=plength(M);
558: for (I=0; I<L; I++) {
559: for (J=0; J<M[I]; J++) {
560: II = I+1; JJ = J+1;
561: V /= (double)bb(M,K,M,II,JJ);
562: /* printf("[%d,%d,%lf]\n",I,J,V);*/
563: }
564: }
565:
566: return(V);
567: }
1.10 takayama 568: static int printp(int kappa[]) {
1.1 takayama 569: int i;
570: printf("(");
571: for (i=0; i<M_n; i++) {
572: if (i <M_n-1) printf("%d,",kappa[i]);
573: else printf("%d)",kappa[i]);
574: }
1.13 takayama 575: return(0);
1.1 takayama 576: }
1.10 takayama 577: static int printp2(int kappa[]) {
1.1 takayama 578: int i,ell;
579: printf("(");
580: ell = plength_t(kappa);
581: for (i=0; i<ell+1; i++) {
582: if (i <ell+1-1) printf("%d,",kappa[i]);
583: else printf("%d)",kappa[i]);
584: }
1.13 takayama 585: return(0);
1.1 takayama 586: }
587:
1.10 takayama 588: static int test_beta() {
1.1 takayama 589: int kappa[M_n0]={2,1,0};
590: int mu1[M_n0]={1,0,0};
591: int mu2[M_n0]={1,1,0};
592: int mu3[M_n0]={2,0,0};
593: printp(kappa); printf(","); printp(mu3); printf(": beta = %lf\n",beta(kappa,mu3));
594: printp(kappa); printf(","); printp(mu1); printf(": beta = %lf\n",beta(kappa,mu1));
595: printp(kappa); printf(","); printp(mu2); printf(": beta = %lf\n",beta(kappa,mu2));
596: }
597:
598: /* main() { test_beta(); } */
599:
600:
601: /*
1.12 takayama 602: cf. w1m.rr
1.1 takayama 603: matrix hypergeometric by jack
604: N variables, up to degree M.
605: */
606: /* todo
1.12 takayama 607: def mhgj(A,B,N,M) {
608: F = 0;
609: P = partition_a(N,M);
610: for (I=0; I<length(P); I++) {
611: K = P[I];
612: F += qk(K,A,B)*jack(K,N);
613: }
614: return(F);
615: }
1.1 takayama 616: */
617:
618:
619: /* The quotient of (3.10) of Koev-Edelman K=kappa, M=mu, SK=k */
1.2 takayama 620: static double q3_10(int K[],int M[],int SK) {
1.1 takayama 621: extern int Alpha;
622: int Mp[M_m_MAX];
1.3 takayama 623: int ML[M_nmx];
624: int N[M_nmx];
1.1 takayama 625: int i,R;
626: double T,Q,V,Ur,Vr,Wr;
627: ptrans(M,Mp);
628: for (i=0; i<M_n; i++) {ML[i] = M[i]; N[i] = M[i];}
629: N[SK-1] = N[SK-1]-1;
630:
631: T = SK-Alpha*M[SK-1];
632: Q = T+1;
633: V = Alpha;
634: for (R=1; R<=SK; R++) {
635: Ur = Q-R+Alpha*K[R-1];
636: V *= Ur/(Ur+Alpha-1);
637: }
638: for (R=1; R<=SK-1; R++) {
639: Vr = T-R+Alpha*M[R-1];
640: V *= (Vr+Alpha)/Vr;
641: }
642: for (R=1; R<=M[SK-1]-1; R++) {
643: Wr = Mp[R-1]-T-Alpha*R;
644: V *= (Wr+Alpha)/Wr;
645: }
646: return(V);
647: }
648:
1.2 takayama 649: static double q3_5(double A[],double B[],int K[],int I) {
1.1 takayama 650: extern int Alpha;
651: int Kp[M_m_MAX];
652: double C,D,V,Ej,Fj,Gj,Hj,Lj;
653: int J,P,Q;
654: ptrans(K,Kp);
655: P=A_LEN;; Q = B_LEN;
656: C = -((double)(I-1))/Alpha+K[I-1]-1;
657: D = K[I-1]*Alpha-I;
658:
659: V=1;
660:
661: for (J=1; J<=P; J++) {
1.12 takayama 662: V *= (A[J-1]+C);
1.1 takayama 663: }
664: for (J=1; J<=Q; J++) {
1.12 takayama 665: V /= (B[J-1]+C);
1.1 takayama 666: }
667:
668: for (J=1; J<=K[I-1]-1; J++) {
1.12 takayama 669: Ej = D-J*Alpha+Kp[J-1];
670: Gj = Ej+1;
671: V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha));
1.1 takayama 672: }
673: for (J=1; J<=I-1; J++) {
1.12 takayama 674: Fj=K[J-1]*Alpha-J-D;
675: Hj=Fj+Alpha;
676: Lj=Hj*Fj;
677: V *= (Lj-Fj)/(Lj+Hj);
1.1 takayama 678: }
679: return(V);
680: }
681:
1.2 takayama 682: static void mtest4() {
1.1 takayama 683: double A[A_LEN] = {1.5};
684: double B[B_LEN]={6.5};
685: int K[M_n0] = {3,2,0};
686: int I=2;
687: int Ki[M_n0]={3,1,0};
688: double V1,V2;
689: V1=q3_5(A,B,K,I);
690: V2=qk(K,A,B)/qk(Ki,A,B);
691: printf("%lf== %lf?\n",V1,V2);
692: }
1.2 takayama 693: static void mtest4b() {
1.1 takayama 694: int K[M_n0]={3,2,0};
695: int M[M_n0]={2,1,0};
696: int N[M_n0]={2,0};
697: int SK=2;
698: double V1,V2;
699: V1=q3_10(K,M,SK);
700: V2=beta(K,N)/beta(K,M);
701: printf("%lf== %lf?\n",V1,V2);
702: }
703:
704: /* main() { mtest4(); mtest4b(); } */
705:
706: /* nk in (4.1),
1.12 takayama 707: */
1.2 takayama 708: static int nk(int KK[]) {
1.1 takayama 709: extern int *Darray;
710: int N,I,Ki;
1.3 takayama 711: int Kpp[M_nmx];
1.1 takayama 712: int i;
713: N = plength(KK);
714: if (N == 0) return(0);
715: if (N == 1) return(KK[0]);
716: for (i=0; i<M_n; i++) Kpp[i] = 0;
717: for (I=0; I<N-1; I++) Kpp[I] = KK[I];
718: Ki = KK[N-1];
719: /* K = (Kpp,Ki) */
720: return(Darray[nk(Kpp)]+Ki-1);
721: }
1.2 takayama 722: static int plength2(int P1[],int P2[]) {
1.1 takayama 723: int S1,S2;
724: S1 = plength(P1); S2 = plength(P2);
725: if (S1 > S2) return(1);
726: else if (S1 == S2) {
727: S1=mysum(P1); S2=mysum(P2);
728: if(S1 > S2) return(1);
729: else if (S1 == S2) return(0);
730: else return(-1);
731: }
732: else return(-1);
733: }
1.2 takayama 734: static int myeq(int P1[],int P2[]) {
1.1 takayama 735: int I,L1;
736: if ((L1=plength(P1)) != plength(P2)) return(0);
737: for (I=0; I<L1; I++) {
738: if (P1[I] != P2[I]) return(0);
739: }
740: return(1);
741: }
742: /*
743: M is a degree, N is a number of variables
744: genDarray(3,3);
745: N(0)=0;
746: N(1)=1;
747: N(2)=2;
748: N(3)=3;
749: N(1,1)=4; D[1] = 4
750: N(2,1)=5; D[2] = 5;
751: N(1,1,1)=6; D[4] = 6;
752: still buggy.
753: */
754:
1.2 takayama 755: static int pListPartition(int M,int N) {
1.1 takayama 756: extern int M_m;
757: extern int M_kap[];
758: int I;
759: /* initialize */
760: if (M_n != N) {
1.3 takayama 761: fprintf(stderr,"M_n != N\n"); mh_exit(-1);
1.1 takayama 762: }
763: M_m = M;
764: /* M_plist = []; */
765: /* end of initialize */
766: (*M_pExec)(); /* exec for 0 */
767: for (I=1; I<=M_n; I++) {
768: pListPartition2(M_m,1,I,M_m);
769: }
770: /* M_plist = reverse(M_plist); */
771: return(1);
772: }
773:
774: /*
775: Enumerate all such that
776: Less >= M_kap[From], ..., M_kap[To], |(M_kap[From],...,M_kap[To])|<=M,
777: */
1.2 takayama 778: static int pListPartition2(int Less,int From,int To, int M) {
1.1 takayama 779: int I,R;
1.11 takayama 780: mh_check_intr(100);
1.1 takayama 781: if (To < From) {
1.12 takayama 782: (*M_pExec)(); return(0);
1.1 takayama 783: }
784: for (I=1; (I<=Less) && (I<=M) ; I++) {
785: M_kap[From-1] = I;
786: R = pListPartition2(I,From+1,To,M-I);
787: }
788: return(1);
789: }
790:
791: /*
792: partition に対してやる仕事をこちらへ書く.
793: */
1.2 takayama 794: static void pExec_0() {
1.1 takayama 795: if (Debug) {
1.12 takayama 796: printf("M_kap=");
797: printp(M_kap);
798: printf("\n");
1.1 takayama 799: }
800: }
801:
802: /* Test.
1.12 takayama 803: Compare pListPartition(4,3); genDarray(4,3);
804: Compare pListPartition(5,3); genDarray(5,3);
1.1 takayama 805:
806: */
807:
808: /*
1.12 takayama 809: main() {
1.1 takayama 810: M_pExec = pExec_0;
811: pListPartition(5,3);
1.12 takayama 812: }
1.1 takayama 813: */
814:
815:
816: /*
817: List all horizontal strips.
818: Kap[0] is not a dummy in C code. !(Start from Kap[1].)
819: */
1.2 takayama 820: static int pListHS(int Kap[],int N) {
1.1 takayama 821: extern int HS_n;
822: extern int HS_mu[];
823: int i;
824: HS_n = N;
825: /* Clear HS_mu. Do not forget when N < M_n */
826: for (i=0; i<M_n; i++) HS_mu[i] = 0;
1.13 takayama 827: return(pListHS2(1,N,Kap));
1.1 takayama 828: }
829:
1.2 takayama 830: static int pListHS2(int From,int To,int Kap[]) {
1.1 takayama 831: int More,I;
832: if (To <From) {(*HS_hsExec)(); return(0);}
833: if (From == HS_n) More=0; else More=Kap[From];
834: for (I=Kap[From-1]; I>= More; I--) {
835: HS_mu[From-1] = I;
836: pListHS2(From+1,To,Kap);
837: }
838: return(1);
839: }
840:
1.2 takayama 841: static void hsExec_0() {
1.1 takayama 842: int i;
843: if(Debug) {printf("hsExec: "); printp(HS_mu); printf("\n");}
844: }
845:
846: /*
847: pListHS([0,4,2,1],3);
848: */
849: /*
1.12 takayama 850: main() {
1.1 takayama 851: int Kap[3]={4,2,1};
852: HS_hsExec = hsExec_0;
853: pListHS(Kap,3);
1.12 takayama 854: }
1.1 takayama 855: */
856:
857: /* The number of partitions <= M, with N parts.
1.12 takayama 858: (0,0,...,0) is excluded.
1.1 takayama 859: */
860: #define aP_pki(i,j) P_pki[(i)*(M+1)+(j)]
1.2 takayama 861: static int pmn(int M,int N) {
1.1 takayama 862: int Min_m_n,I,K,S,T,i,j;
863: extern int P_pmn;
864: extern int *P_pki;
865: Min_m_n = (M>N?N:M);
866: /* P_pki=newmat(Min_m_n+1,M+1); */
867: P_pki = (int *) mymalloc(sizeof(int)*(Min_m_n+1)*(M+1));
868: for (i=0; i<Min_m_n+1; i++) for (j=0; j<M+1; j++) aP_pki(i,j) = 0;
869: for (I=1; I<=M; I++) aP_pki(1,I) = 1;
870: for (K=1; K<=Min_m_n; K++) aP_pki(K,0) = 0;
871: S = M;
872: for (K=2; K<=Min_m_n; K++) {
873: for (I=1; I<=M; I++) {
874: if (I-K < 0) T=0; else T=aP_pki(K,I-K);
875: aP_pki(K,I) = aP_pki(K-1,I-1)+T;
876: S += aP_pki(K,I);
877: }
878: }
879: P_pmn=S;
880: if (Debug) {
881: printf("P_pmn=%d\n",P_pmn);
882: for (i=0; i<=Min_m_n; i++) {
883: for (j=0; j<=M; j++) printf("%d,",aP_pki(i,j));
884: printf("\n");
885: }
886: }
887: myfree(P_pki); P_pki=NULL;
888: return(S);
889: }
890:
891: /*
1.12 takayama 892: main() {pmn(4,3); printf("P_pmn=%d\n",P_pmn);}
1.1 takayama 893: */
894:
1.2 takayama 895: static int *cloneP(int a[]) {
1.1 takayama 896: int *p;
897: int i;
898: p = (int *) mymalloc(sizeof(int)*M_n);
899: for (i=0; i<M_n; i++) p[i] = a[i];
900: return(p);
901: }
1.10 takayama 902: static int copyP(int p[],int a[]) {
1.1 takayama 903: int i;
904: for (i=0; i<M_n; i++) p[i] = a[i];
1.13 takayama 905: return(0);
1.1 takayama 906: }
907:
1.2 takayama 908: static void pExec_darray(void) {
1.1 takayama 909: extern int DR_parray;
910: extern int M_kap[];
911: extern int **Parray;
912: extern int *ParraySize;
913: int *K;
914: pExec_0();
915: K = cloneP(M_kap);
916: Parray[DR_parray] = K;
917: ParraySize[DR_parray] = mysum(K);
918: DR_parray++;
919: }
1.10 takayama 920: static int genDarray2(int M,int N) {
1.1 takayama 921: extern int *Darray;
922: extern int **Parray;
923: extern int DR_parray;
924: extern int M_m;
925: int Pmn,I,J,Ksize,i;
926: int **L;
927: int *Nk;
928: int *K;
1.3 takayama 929: int Kone[M_nmx];
1.1 takayama 930:
931: M_m = M;
932: Pmn = pmn(M,N)+1;
933: if (Debug) printf("Degree M = %d, N of vars N = %d, Pmn+1=%d\n",M,N,Pmn);
934: Darray=(int *) mymalloc(sizeof(int)*Pmn);
935: for (i=0; i<Pmn; i++) Darray[i] = 0;
936: Parray=(int **) mymalloc(sizeof(int *)*Pmn);
937: for (i=0; i<Pmn; i++) Parray[i] = NULL;
938: ParraySize=(int *) mymalloc(sizeof(int *)*Pmn);
939: for (i=0; i<Pmn; i++) ParraySize[i] = 0;
940: DR_parray=0;
941: M_pExec = pExec_darray;
942: pListPartition(M,N); /* pExec_darray() is executed for all partitions */
943: L = Parray;
944:
945: Nk = (int *) mymalloc(sizeof(int)*(Pmn+1));
946: for (I=0; I<Pmn; I++) Nk[I] = I;
947: for (I=0; I<Pmn; I++) {
1.11 takayama 948: mh_check_intr(100);
1.12 takayama 949: K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */
950: Ksize = plength(K);
951: if (Ksize >= M_n) {
952: if (Debug) {fprintf(stderr,"Ksize >= M_n\n");}
953: continue;
954: }
955: for (i=0; i<M_n; i++) Kone[i] = 0;
956: for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1;
957: for (J=0; J<Pmn; J++) {
958: if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */
959: }
1.1 takayama 960: }
961: if (Debug) {
1.12 takayama 962: printf("Darray=\n");
963: for (i=0; i<Pmn; i++) printf("%d\n",Darray[i]);
964: printf("-----------\n");
1.1 takayama 965: }
1.13 takayama 966: return(0);
1.1 takayama 967: }
968:
969:
970: /* main() { genDarray2(4,3);} */
971:
972: /* M_beta_0[*] value of beta_{kappa,mu}, M_beta_1[*] N_mu */
1.10 takayama 973: static int isHStrip(int Kap[],int Nu[]) {
1.1 takayama 974: int N1,N2,I,P;
975: N1 = plength(Kap); N2 = plength(Nu);
976: if (N2 > N1) return(0);
977: for (I=0; I<N2; I++) {
978: if (I >= N1-1) P = 0; else P=Kap[I+1];
979: if (Kap[I] < Nu[I]) return(0);
980: if (Nu[I] < P) return(0);
981: }
982: return(1);
983: }
984:
1.2 takayama 985: static void hsExec_beta(void) {
1.1 takayama 986: int *Mu;
987: int N,Nmu,Nnu,Done,J,K,OK,I,RR;
988: int Kapt[M_m_MAX];
989: int Nut[M_m_MAX];
1.3 takayama 990: int Nu[M_nmx];
1.1 takayama 991: int rrMax;
992: hsExec_0();
993: /* printf("M_beta_pt=%a\n",M_beta_pt); */
994: /* Mu = cdr(vtol(HS_mu)); */
995: Mu = HS_mu; /* buggy? need cloneP */
996: if (M_beta_pt == 0) {
997: M_beta_0[0] = 1; M_beta_1[0] = nk(Mu);
998: M_beta_pt++; return;
999: }
1000:
1001: N = HS_n;
1002: Nmu = nk(Mu);
1003: M_beta_1[M_beta_pt] = Nmu;
1004: ptrans(M_beta_kap,Kapt);
1005: /* Mu, Nu is exchanged in this code. cf. the K-E paper */
1006: copyP(Nu,Mu); /* buggy need clone? */
1007: for (I=0; I<N; I++) {
1008: Nu[I]++;
1009: if (!isHStrip(M_beta_kap,Nu)) {Nu[I]--; continue;}
1010: Nnu = nk(Nu);
1011: ptrans(Nu,Nut);
1012: Done=0;
1013: for (J=M_beta_pt-1; J>=0; J--) {
1014: if (M_beta_1[J] == Nnu) {
1.12 takayama 1015: K=I+1;
1016: if (Debug) {
1017: printf("Found at J=%d, K=%d, q3_10(Kap,Nu,K)=%lf,Nu,Mu= \n",
1018: J,K,q3_10(M_beta_kap,Nu,K));
1019: printp(Nu); printf("\n");
1020: printp(Mu); printf("\n");
1021: }
1022: /* Check other conditions. See Numata's mail on Dec 24, 2011. */
1023: rrMax = Nu[I]-1;
1024: if ((plength_t(Kapt) < rrMax) || (plength_t(Nut) < rrMax)) {
1025: if (Debug) printf(" is not taken (length). \n");
1026: break;
1027: }
1028: OK=1;
1029: for (RR=0; RR<rrMax; RR++) {
1030: if (Kapt[RR] != Nut[RR]) { OK=0; break;}
1031: }
1032: if (!OK) { if (Debug) printf(" is not taken.\n"); break; }
1033: /* check done. */
1034: M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K);
1035: Done = 1; break;
1.1 takayama 1036: }
1037: }
1038: if (Done) break; else Nu[I]--;
1039: }
1040: if (!Done) {
1041: if (Debug) printf("BUG: not found M_beta_pt=%d.\n",M_beta_pt);
1.13 takayama 1042: /* M_beta_0[M_beta_pt] = NAN; error("Not found."); */
1.1 takayama 1043: M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu);
1044: }
1045: /* Fix the bug of mh.rr */
1046: M_beta_pt++;
1047: }
1.10 takayama 1048: static int genBeta(int Kap[]) {
1.1 takayama 1049: extern double *M_beta_0;
1050: extern int *M_beta_1;
1051: extern int M_beta_pt;
1052: extern int M_beta_kap[];
1053: extern int P_pmn;
1054: int I,J,N;
1055: if (Debug) {printp(Kap); printf("<-Kappa, P_pmn=%d\n",P_pmn);}
1056: /* M_beta = newmat(2,P_pmn+1); */
1057: M_beta_0 = (double *)mymalloc(sizeof(double)*(P_pmn+1));
1.8 takayama 1058: M_beta_1 = (int *)mymalloc(sizeof(int)*(P_pmn+1));
1.1 takayama 1059: M_beta_pt = 0;
1060: for (I=0; I<=P_pmn; I++) {M_beta_0[I] = NAN; M_beta_1[I] = -1;}
1061: N = plength(Kap);
1062: HS_hsExec = hsExec_beta;
1063: copyP(M_beta_kap,Kap);
1064: pListHS(Kap,N);
1065: }
1066: /*
1067: genDarray2(4,3);
1068: genBeta([2,2,0]);
1069: genBeta([2,1,1]);
1070: */
1071:
1.10 takayama 1072: static int checkBeta1() {
1.1 takayama 1073: int Kap[3] = {2,2,0};
1074: int Kap2[3] = {2,1,0};
1075: int I;
1076: int *Mu;
1077: double Beta_km;
1078: genDarray2(4,3);
1079: genBeta(Kap);
1080: for (I=0; I<M_beta_pt; I++) {
1081: Mu = Parray[M_beta_1[I]];
1082: Beta_km = M_beta_0[I];
1.12 takayama 1083: if (Debug) {
1084: printp(Kap); printf("<--Kap, ");
1085: printp(Mu); printf("<--Mu,");
1086: printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap,Mu));
1087: }
1.1 takayama 1088: }
1089: if (Debug) printf("-------------------------------------\n");
1090: genBeta(Kap2);
1091: for (I=0; I<M_beta_pt; I++) {
1092: Mu = Parray[M_beta_1[I]];
1093: Beta_km = M_beta_0[I];
1.12 takayama 1094: if (Debug) {
1095: printp(Kap2); printf("<--Kap, ");
1096: printp(Mu); printf("<--Mu,");
1097: printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap2,Mu));
1098: }
1.1 takayama 1099: }
1.13 takayama 1100: return(0);
1.1 takayama 1101: }
1102:
1103: /*
1.12 takayama 1104: def checkBeta2() {
1.1 takayama 1105: genDarray2(3,3);
1106: Kap = [2,1,0];
1107: printf("Kap=%a\n",Kap);
1108: genBeta(Kap);
1109: for (I=0; I<M_beta_pt; I++) {
1.12 takayama 1110: Mu = Parray[M_beta[1][I]];
1111: Beta_km = M_beta[0][I];
1112: printf("Mu=%a,",Mu);
1113: printf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));
1114: }
1.1 takayama 1115: }
1116: */
1117:
1118: /* main() { checkBeta1(); } */
1119:
1.2 takayama 1120: static int psublen(int Kap[],int Mu[]) {
1.1 takayama 1121: int L1,L2,A,I;
1122: L1 = plength(Kap);
1123: L2 = plength(Mu);
1124: if (L2 > L1) myerror("psub, length mismatches.");
1125: A = 0;
1126: for (I=0; I<L2; I++) {
1127: if (Kap[I] < Mu[I]) myerror("psub, not Kap >= Mu");
1128: A += Kap[I]-Mu[I];
1129: }
1130: for (I=L2; I<L1; I++) A += Kap[I];
1131: return(A);
1132: }
1133:
1134:
1135: /* Table of Jack polynomials
1.12 takayama 1136: Jack[1]* one variable.
1137: Jack[2]* two variables.
1.1 takayama 1138: ...
1.12 takayama 1139: Jack[M_n]* n variables.
1140: Jack[P][J]*
1141: D^J(P variables jack of p variables). Example. J=001 d_1, 010 d_2, 100 d_3
1142: 0<=J<=2^{M_n}-1
1143: Jack[P][J][nk(Kappa)] Jack_Kappa, Kappa is a partition.
1144: 0<=nk(Kappa)<=pmn(M_m,M_n)
1.1 takayama 1145: */
1146:
1147: #define aM_jack(i,j,k) ((M_jack[i])[(j)*(Pmn+1)+(k)])
1.10 takayama 1148: static int genJack(int M,int N) {
1.1 takayama 1149: extern double **M_jack;
1150: extern int M_2n;
1151: extern int P_pmn;
1152: extern int *M_beta_1;
1153: int Pmn,I,J,K,L,Nv,H,P;
1154: int *Kap,*Mu;
1155: double Jack,Beta_km;
1156: int Nk,JJ;
1.5 takayama 1157: if (Debug) printf("genJack(%d,%d)\n",M,N);
1.3 takayama 1158: M_jack = (double **) mymalloc(sizeof(double *)*(N+2));
1.1 takayama 1159: M_2n = imypower(2,N);
1160: Pmn = pmn(M,N); /*P_pmn is initializeded.
1161: Warning. It is reset when pmn is called.*/
1162: for (I=0; I<=N; I++) M_jack[I] = (double *)mymalloc(sizeof(double)*(M_2n*(Pmn+1))); /* newmat(M_2n,Pmn+1); */
1.3 takayama 1163: M_jack[N+1] = NULL;
1.1 takayama 1164: genDarray2(M,N); /* Darray, Parray is initialized */
1165: for (I=1; I<=N; I++) aM_jack(I,0,0) = 1;
1166: if (M_df) {
1167: for (I=1; I<=N; I++) {
1168: for (J=1; J<M_2n; J++) aM_jack(I,J,0) = 0;
1169: }
1170: }
1171:
1.3 takayama 1172: /* N must satisfies N > 0 */
1.1 takayama 1173: for (K=1; K<=M; K++) {
1174: aM_jack(1,0,K) = jack1(K);
1175: if (M_df) {
1176: aM_jack(1,1,K) = jack1diff(K); /* diff(jack([K],1),x_1); */
1177: for (J=2; J<M_2n; J++) aM_jack(1,J,K) = 0;
1178: }
1179: }
1180: for (I=1; I<=N; I++) {
1181: for (K=M+1; K<Pmn+1; K++) {
1182: aM_jack(I,0,K) = NAN;
1183: if (M_df) {
1184: for (J=1; J<M_2n; J++) {
1185: if (J >= 2^I) aM_jack(I,J,K) = 0;
1.12 takayama 1186: else aM_jack(I,J,K) = NAN;
1.1 takayama 1187: }
1188: }
1189: }
1190: }
1191:
1192: /* Start to evaluate the entries of the table */
1193: for (K=1; K<=Pmn; K++) {
1194: Kap = Parray[K]; /* bug. need copy? */
1195: L = plength(Kap);
1196: for (I=1; I<=L-1; I++) {
1197: aM_jack(I,0,K) = 0;
1198: if (M_df) {
1199: for (J=1; J<M_2n; J++) aM_jack(I,J,K) = 0;
1200: }
1201: }
1202: if (Debug) {printf("Kappa="); printp(Kap);}
1203: /* Enumerate horizontal strip of Kappa */
1204: genBeta(Kap); /* M_beta_pt stores the number of hs */
1205: /* Nv is the number of variables */
1206: for (Nv = (L==1?2:L); Nv <= N; Nv++) {
1207: Jack = 0;
1208: for (H=0; H<M_beta_pt; H++) {
1209: Nk = M_beta_1[H];
1210: Mu = Parray[Nk];
1.12 takayama 1211: if (UseTable) {
1212: Beta_km = M_beta_0[H];
1213: }else{
1214: Beta_km = beta(Kap,Mu);
1215: /* do not use the M_beta table. It's buggy. UseTable is experimental.*/
1216: }
1.1 takayama 1217: if (Debug) {printf("Nv(number of variables)=%d, Beta_km=%lf, Mu=",Nv,Beta_km);
1.12 takayama 1218: printp(Mu); printf("\n");}
1.1 takayama 1219: P = psublen(Kap,Mu);
1220: Jack += aM_jack(Nv-1,0,Nk)*Beta_km*xval(Nv,P); /* util_v(x,[Nv])^P;*/
1.12 takayama 1221: if (Debug) printf("xval(%d,%d)=%lf\n",Nv,P,xval(Nv,P));
1.1 takayama 1222: }
1223: aM_jack(Nv,0,K) = Jack;
1224: if (M_df) {
1225: /* The case of M_df > 0. */
1226: for (J=1; J<M_2n; J++) {
1.12 takayama 1227: mh_check_intr(100);
1228: Jack = 0;
1229: for (H=0; H<M_beta_pt; H++) {
1230: Nk = M_beta_1[H];
1231: Mu = Parray[Nk];
1232: if (UseTable) {
1233: Beta_km = M_beta_0[H];
1234: }else{
1235: Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */
1236: }
1237: if (Debug) {printf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);
1238: printp(Mu); printf("\n"); }
1239: P = psublen(Kap,Mu);
1240: if (J & (1 << (Nv-1))) {
1241: JJ = J & ((1 << (Nv-1)) ^ 0xffff); /* NOTE!! Up to 16 bits. mh-15 */
1242: if (P != 0) {
1243: Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1);
1244: }
1245: }else{
1246: Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P);
1247: }
1248: }
1249: aM_jack(Nv,J,K) = Jack;
1250: if (Debug) printf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);
1.1 takayama 1251: } /* end of J loop */
1252: }
1253: }
1254: }
1255: }
1256:
1257:
1258: /* checkJack1(3,3)
1.12 takayama 1259: */
1.10 takayama 1260: static int checkJack1(int M,int N) {
1.1 takayama 1261: int I,K;
1262: extern int P_pmn;
1263: extern double M_x[];
1264: int Pmn; /* used in aM_jack */
1265: /* initialize x vars. */
1266: for (I=1; I<=N; I++) {
1267: M_x[I-1] = ((double)I)/10.0;
1268: }
1269: genJack(M,N);
1270: Pmn = P_pmn;
1271: for (I=1; I<=N; I++) {
1272: for (K=0; K<=P_pmn; K++) {
1273: printp(Parray[K]);
1274: printf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
1275: }
1276: }
1277: for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);
1278: printf("<--x\n");
1.13 takayama 1279: return(0);
1.1 takayama 1280: }
1281: /*main() { checkJack1(3,3); }*/
1282:
1283:
1.10 takayama 1284: static int checkJack2(int M,int N) {
1.1 takayama 1285: int I,K,J;
1286: extern int P_pmn;
1287: extern double M_x[];
1.10 takayama 1288: extern int M_df;
1.1 takayama 1289: int Pmn; /* used in aM_jack */
1290: M_df=1;
1291: /* initialize x vars. */
1292: for (I=1; I<=N; I++) {
1293: M_x[I-1] = ((double)I)/10.0;
1294: }
1295: genJack(M,N);
1296: Pmn = P_pmn;
1297: for (I=1; I<=N; I++) {
1298: for (K=0; K<=P_pmn; K++) {
1299: printp(Parray[K]);
1300: printf("<--Kap, Nv=%d, TableJack=%lf\n",I,aM_jack(I,0,K));
1301: }
1302: }
1303: for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);
1304: printf("<--x\n");
1305:
1306: for (I=1; I<=N; I++) {
1307: for (K=0; K<=P_pmn; K++) {
1308: for (J=0; J<M_2n; J++) {
1.12 takayama 1309: printp(Parray[K]);
1310: printf("<--Kap, Nv=%d,J(diff)=%d, D^J Jack=%lf\n",
1311: I,J,aM_jack(I,J,K));
1312: }
1313: }
1.1 takayama 1314: }
1.13 takayama 1315: return(0);
1.1 takayama 1316: }
1317:
1318: /* main() { checkJack2(3,3); } */
1319:
1320: double mh_t(double A[A_LEN],double B[B_LEN],int N,int M) {
1321: double F,F2;
1322: extern int M_df;
1323: extern int P_pmn;
1324: extern double *M_qk;
1325: extern double M_rel_error;
1.15 takayama 1326: extern int M_m;
1.16 takayama 1327: extern int M_m_estimated_approx_deg;
1328: extern double M_assigned_series_error;
1.1 takayama 1329: int Pmn;
1330: int K;
1331: int *Kap;
1332: int size;
1.15 takayama 1333: int i;
1334: double partial_sum[M_m_MAX+1];
1.16 takayama 1335: double iv;
1336: double serror;
1.1 takayama 1337: F = 0; F2=0;
1338: M_df=1;
1339: genJack(M,N);
1.8 takayama 1340: M_qk = (double *)mymalloc(sizeof(double)*(P_pmn+1)); /* found a bug. */
1.1 takayama 1341: Pmn = P_pmn;
1342: size = ParraySize[P_pmn];
1343: for (K=0; K<=P_pmn; K++) {
1.11 takayama 1344: mh_check_intr(100);
1.12 takayama 1345: Kap = Parray[K];
1346: M_qk[K] = qk(Kap,A,B);
1347: F += M_qk[K]*aM_jack(N,0,K);
1348: if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K);
1349: if (Debug) printf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size);
1350: if (Debug && (ParraySize[K] == size)) printf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K));
1.1 takayama 1351: }
1352: M_rel_error = F-F2;
1.15 takayama 1353:
1.16 takayama 1354: M_m_estimated_approx_deg = -1; serror=1;
1.15 takayama 1355: for (i=0; i<=M_m; i++) {
1356: partial_sum[i] = 0.0; partial_sum[i+1] = 0.0;
1357: for (K=0; K<=P_pmn; K++) {
1358: if (ParraySize[K] == i) partial_sum[i] += M_qk[K]*aM_jack(N,0,K);
1359: }
1360: if (i>0) partial_sum[i] += partial_sum[i-1];
1.16 takayama 1361: serror = myabs((partial_sum[i]-partial_sum[i-1])/partial_sum[i-1]);
1362: if ((i>0)&&(M_m_estimated_approx_deg < 0)&&(serror<M_assigned_series_error)) {
1363: M_m_estimated_approx_deg = i; break;
1364: }
1365: }
1366: if (M_m_estimated_approx_deg < 0) {
1367: M_m_estimated_approx_deg = M_m+mymin(5,mymax(1,(int)log(serror/M_assigned_series_error))); /* Heuristic */
1.15 takayama 1368: }
1369: /*
1370: for (K=0; K<=P_pmn; K++) {
1371: printf("Kappa="); for (i=0; i<N; i++) printf("%d ",Parray[K][i]); printf("\n");
1372: printf("ParraySize(%d)=%d (|kappa|), M_m=%d\n",K,ParraySize[K],M_m);
1373: }
1374: for (i=0; i<=M_m; i++) {
1375: printf("partial_sum[%d]=%lg\n",i,partial_sum[i]);
1376: }
1377: */
1.16 takayama 1378: M_estimated_X0g = X0g;
1379: iv=myabs(F*iv_factor());
1380: if (iv < M_x0value_min) M_estimated_X0g = X0g*mymax(2,log(log(1/iv))); /* This is heuristic */
1.19 takayama 1381: M_estimated_X0g = mymin(M_estimated_X0g,M_X0g_bound);
1.16 takayama 1382: M_mh_t_success = 1;
1.19 takayama 1383: if (M_estimated_X0g != X0g) M_mh_t_success=0;
1384: if (M_m_estimated_approx_deg > M_m) M_mh_t_success=0;
1.16 takayama 1385:
1.19 takayama 1386: M_series_error = serror;
1387: M_recommended_abserr = iv*M_assigned_series_error;
1.27 takayama 1388: M_recommended_relerr = M_series_error;
1.25 takayama 1389:
1390: if (M_show_autosteps) {
1391: printf("%%%%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);
1392: printf("%%%%x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
1393: printf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
1394: fprintf(stderr,"%%%%(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);
1395: fprintf(stderr,"%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
1396: fprintf(stderr,"%%%%(stderr) F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
1397: }
1.16 takayama 1398:
1.19 takayama 1399: M_mh_t_value=F;
1.1 takayama 1400: return(F);
1401: }
1402: double mh_t2(int J) {
1403: extern double *M_qk;
1404: double F;
1405: int K;
1406: int Pmn;
1407: extern int P_pmn;
1.3 takayama 1408: if (M_qk == NULL) {myerror("Call mh_t first."); mh_exit(-1); }
1.1 takayama 1409: F = 0;
1410: Pmn = P_pmn;
1411: for (K=0; K<P_pmn; K++) {
1.12 takayama 1412: F += M_qk[K]*aM_jack(M_n,J,K);
1.1 takayama 1413: }
1414: return(F);
1415: }
1416:
1.10 takayama 1417: static int mtest1b() {
1.1 takayama 1418: double A[1] = {1.5};
1419: double B[1] = {1.5+5};
1420: int I,N,M,J;
1421: double F;
1422: N=3; M=6;
1423: for (I=1; I<=N; I++) {
1424: M_x[I-1] = ((double)I)/10.0;
1425: }
1426: mh_t(A,B,N,M);
1427: for (J=0; J<M_2n; J++) {
1.12 takayama 1428: F=mh_t2(J);
1429: printf("J=%d, D^J mh_t=%lf\n",J,F);
1.1 takayama 1430: }
1431: }
1432:
1433: /* main() { mtest1b(); }*/
1434:
1435:
1436:
1437:
1438: #define TEST 1
1439: #ifndef TEST
1440:
1441: #endif
1442:
1443: /****** from mh-n.c *****/
1444:
1445: #define SMAX 4096
1.3 takayama 1446: #define inci(i) { i++; if (i >= argc) { fprintf(stderr,"Option argument is not given.\n"); return(NULL); }}
1.1 takayama 1447:
1.2 takayama 1448: static int imypower(int x,int n) {
1.1 takayama 1449: int i;
1450: int v;
1.3 takayama 1451: if (n < 0) {myerror("imypower"); mh_exit(-1);}
1.1 takayama 1452: v = 1;
1453: for (i=0; i<n; i++) v *= x;
1454: return(v);
1455: }
1456:
1.11 takayama 1457: #ifdef STANDALONE2
1.1 takayama 1458: main(int argc,char *argv[]) {
1.3 takayama 1459: mh_exit(MH_RESET_EXIT);
1.12 takayama 1460: /* jk_main(argc,argv);
1461: printf("second run.\n"); */
1.2 takayama 1462: jk_main(argc,argv);
1.17 takayama 1463: return(0);
1.2 takayama 1464: }
1.3 takayama 1465: #endif
1466:
1467: struct MH_RESULT *jk_main(int argc,char *argv[]) {
1.16 takayama 1468: int i;
1469: struct MH_RESULT *ans;
1470: extern int M_automatic;
1471: extern int M_mh_t_success;
1472: extern double M_estimated_X0g;
1473: extern int M_m_estimated_approx_deg;
1474: for (i=1; i<argc; i++) {
1475: if (strcmp(argv[i],"--automatic")==0) {
1476: inci(i);
1477: sscanf(argv[i],"%d",&M_automatic);
1478: break;
1479: }
1480: }
1481: ans=jk_main2(argc,argv,0,0.0,0);
1482: if (!M_automatic) return(ans);
1483: if (M_mh_t_success) return(ans);
1484: while (!M_mh_t_success) {
1485: ans=jk_main2(argc,argv,1,M_estimated_X0g,M_m_estimated_approx_deg);
1486: }
1487: return(ans);
1488: }
1489:
1490: struct MH_RESULT *jk_main2(int argc,char *argv[],int automode,double newX0g,int newDegree) {
1.1 takayama 1491: double *y0;
1492: double x0,xn;
1493: double ef;
1494: int i,j,rank;
1495: double a[1]; double b[1];
1496: extern double M_x[];
1497: extern double *Beta;
1498: extern int M_2n;
1.16 takayama 1499: extern int M_mh_t_success;
1.3 takayama 1500: char swork[1024];
1501: struct MH_RESULT *ans=NULL;
1502: struct SFILE *ofp = NULL;
1503: int idata=0;
1504: JK_byFile = 1;
1505: jk_initializeWorkArea();
1.1 takayama 1506: UseTable = 1;
1.3 takayama 1507: Mapprox=6;
1.1 takayama 1508: for (i=1; i<argc; i++) {
1.12 takayama 1509: if (strcmp(argv[i],"--idata")==0) {
1510: inci(i);
1511: setParam(argv[i]); idata=1;
1512: }else if (strcmp(argv[i],"--degree")==0) {
1513: inci(i);
1514: sscanf(argv[i],"%d",&Mapprox);
1515: }else if (strcmp(argv[i],"--x0")==0) {
1516: inci(i);
1517: sscanf(argv[i],"%lg",&X0g);
1518: }else if (strcmp(argv[i],"--notable")==0) {
1519: UseTable = 0;
1520: }else if (strcmp(argv[i],"--debug")==0) {
1521: Debug = 1;
1522: }else if (strcmp(argv[i],"--help")==0) {
1523: usage(); return(0);
1524: }else if (strcmp(argv[i],"--bystring")==0) {
1525: if (idata) {fprintf(stderr,"--bystring must come before --idata option.\n"); mh_exit(-1);}
1526: JK_byFile = 0;
1.16 takayama 1527: }else if (strcmp(argv[i],"--automatic")==0) {
1528: inci(i); /* ignore, in this function */
1.23 takayama 1529: }else if (strcmp(argv[i],"--assigned_series_error")==0) {
1.16 takayama 1530: inci(i);
1531: sscanf(argv[i],"%lg",&M_assigned_series_error);
1532: }else if (strcmp(argv[i],"--x0value_min")==0) {
1533: inci(i);
1534: sscanf(argv[i],"%lg",&M_x0value_min);
1.12 takayama 1535: }else {
1536: fprintf(stderr,"Unknown option %s\n",argv[i]);
1537: usage();
1538: return(NULL);
1539: }
1.1 takayama 1540: }
1.16 takayama 1541: if (!idata) setParam(NULL);
1542: if (automode) {
1543: Mapprox = newDegree;
1544: X0g = newX0g;
1545: }
1.3 takayama 1546:
1547: /* Initialize global variables */
1548: M_n = Mg;
1549: HS_n=M_n;
1.16 takayama 1550: if (!JK_byFile) {
1551: ans = (struct MH_RESULT *)mymalloc(sizeof(struct MH_RESULT));
1552: ans->message = NULL;
1.19 takayama 1553: ans->t_success = 0;
1554: ans->series_error = 1.0e+10;
1555: ans->recommended_abserr = 1.0e-10;
1.16 takayama 1556: }
1.3 takayama 1557: else ans = NULL;
1.26 takayama 1558: if (M_automatic) {
1559: /* Differentiation can be M_m in the bit pattern in the M_n variable case.*/
1560: if (M_n > Mapprox) Mapprox=M_n;
1561: }
1.3 takayama 1562: /* Output by a file=stdout */
1563: ofp = mh_fopen("stdout","w",JK_byFile);
1.1 takayama 1564:
1.3 takayama 1565: sprintf(swork,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp);
1566: sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp);
1.1 takayama 1567: if (M_n != Mg) {
1.12 takayama 1568: myerror("Mg must be equal to M_n\n"); mh_exit(-1);
1.1 takayama 1569: }
1.8 takayama 1570: if (Debug) showParam(NULL,1);
1.1 takayama 1571: for (i=0; i<M_n; i++) {
1.12 takayama 1572: M_x[i] = Beta[i]*X0g;
1.1 takayama 1573: }
1.19 takayama 1574:
1575: M_beta_i_x_o2_max=myabs(M_x[0]/2);
1576: if (M_n <= 1) M_beta_i_beta_j_min = myabs(Beta[0]);
1577: else M_beta_i_beta_j_min = myabs(Beta[1]-Beta[0]);
1578: for (i=0; i<M_n; i++) {
1579: if (myabs(M_x[i]/2) > M_beta_i_x_o2_max) M_beta_i_x_o2_max = myabs(M_x[i]/2);
1580: for (j=i+1; j<M_n; j++) {
1581: if (myabs(Beta[i]-Beta[j]) < M_beta_i_beta_j_min)
1582: M_beta_i_beta_j_min = myabs(Beta[i]-Beta[j]);
1583: }
1584: }
1585:
1.1 takayama 1586: a[0] = ((double)Mg+1.0)/2.0;
1587: b[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */
1588: if (Debug) printf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",a[0],b[0],M_n,Mapprox);
1589: mh_t(a,b,M_n,Mapprox);
1.19 takayama 1590: if ((!M_mh_t_success) && M_automatic) {
1.16 takayama 1591: jk_freeWorkArea();
1592: return NULL;
1593: }
1.1 takayama 1594: if (imypower(2,M_n) != M_2n) {
1.12 takayama 1595: sprintf(swork,"M_n=%d,M_2n=%d\n",M_n,M_2n); mh_fputs(swork,ofp);
1596: myerror("2^M_n != M_2n\n"); mh_exit(-1);
1.1 takayama 1597: }
1.3 takayama 1598: sprintf(swork,"%%%%M_rel_error=%lg\n",M_rel_error); mh_fputs(swork,ofp);
1.1 takayama 1599: for (j=0; j<M_2n; j++) {
1.12 takayama 1600: Iv[j] = mh_t2(j);
1.1 takayama 1601: }
1602: Ef = iv_factor();
1.8 takayama 1603: showParam(ofp,0);
1.3 takayama 1604:
1605: /* return the result */
1.16 takayama 1606: if (!JK_byFile) {
1.12 takayama 1607: ans->x = X0g;
1608: ans->rank = imypower(2,Mg);
1609: ans->y = (double *)mymalloc(sizeof(double)*(ans->rank));
1610: for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i];
1611: ans->size=1;
1612: ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size));
1613: (ans->sfpp)[0] = ofp;
1.19 takayama 1614:
1615: ans->t_success = M_mh_t_success;
1616: ans->series_error = M_series_error;
1617: ans->recommended_abserr = M_recommended_abserr;
1.3 takayama 1618: }
1.7 takayama 1619: if (Debug) printf("jk_freeWorkArea() starts\n");
1.3 takayama 1620: jk_freeWorkArea();
1.7 takayama 1621: if (Debug) printf("jk_freeWorkArea() has finished.\n");
1.3 takayama 1622: return(ans);
1.1 takayama 1623: }
1624:
1.10 takayama 1625: static int usage() {
1.1 takayama 1626: fprintf(stderr,"Usages:\n");
1.14 takayama 1627: fprintf(stderr,"hgm_jack-n [--idata input_data_file --x0 x0 --degree approxm]\n");
1.23 takayama 1628: fprintf(stderr," [--automatic n --assigned_series_error e --x0value_min e2]\n");
1.28 ! takayama 1629: fprintf(stderr,"\nThe command hgm_jack-n [options] generates an input for hgm_w-n, Pr({y | y<xmax}), which is the cumulative distribution function of the largest root of the m by m Wishart matrices with n degrees of freedom and the covariantce matrix sigma.\n");
1.14 takayama 1630: fprintf(stderr,"The hgm_jack-n uses the Koev-Edelman algorithm to evalute the matrix hypergeometric function.\n");
1.1 takayama 1631: fprintf(stderr,"The degree of the approximation (Mapprox) is given by the --degree option.\n");
1.28 ! takayama 1632: fprintf(stderr,"Parameters are specified by the input_data_file. Otherwise, default values are used.\n\n");
! 1633: fprintf(stderr,"The format of the input_data_file: (The orders of the input data must be kept.)\n");
1.1 takayama 1634: fprintf(stderr," Mg: m(the number of variables), Beta: beta=sigma^(-1)/2 (diagonized), Ng: n,\n");
1.24 takayama 1635: fprintf(stderr," (Add a comment line %%Ng= before the data Ng to check the number of beta.)\n");
1.14 takayama 1636: fprintf(stderr," X0g: starting value of x(when --x0 option is used, this value is used)\n");
1637: fprintf(stderr," Iv: initial values at X0g*Beta (see our paper how to order them), are evaluated in this program. Give zeros or the symbol * to skip rank many inputs.\n");
1.1 takayama 1638: fprintf(stderr," Ef: a scalar factor to the initial value. It is calculated by this program. Give the zero.\n");
1.14 takayama 1639: fprintf(stderr," Hg: h (step size) which is for hgm_w-n, \n");
1.28 ! takayama 1640: fprintf(stderr," Dp: output data is stored in every Dp steps when output_data_file is specified. This is for hgm_w-n.\n");
! 1641: fprintf(stderr," Xng: terminating value of x. This is for hgm_w-n.\n");
! 1642: fprintf(stderr,"Optional parameters automatic, ... are interpreted by a parser. See setParam() in jack-n.c and Testdata/tmp-idata2.txt as an example. Optional paramters are given as %%parameter_name=value Lines starting with %%%% or # are comment lines.\n");
! 1643: fprintf(stderr,"Parameters are redefined when they appear more than once in the idata file and the command line options.\n\n");
! 1644: fprintf(stderr,"With the --notable option, it does not use the Lemma 3.2 of Koev-Edelman (there is a typo: kappa'_r = mu'_r for 1<=r<=mu_k).\n");
! 1645: fprintf(stderr,"An example format of the input_data_file can be obtained by executing hgm_jack-n with no option. When there is no --idata file, all options are ignored.\n");
! 1646: fprintf(stderr,"By --automatic option, X0g and degree are automatically determined from assigend_series_error. The current strategy is described in mh_t in jack-n.c\n");
1.26 takayama 1647: fprintf(stderr,"Default values for the papameters of the automatic mode: automatic=%d, assigned_series_error=%lg, x0value_min=%lg\n",M_automatic,M_assigned_series_error,M_x0value_min);
1.20 takayama 1648: fprintf(stderr,"Todo: automatic mode throws away the table of Jack polynomials of the previous degrees and reevaluate them. They should be kept.\n");
1.1 takayama 1649: fprintf(stderr,"\nExamples:\n");
1.14 takayama 1650: fprintf(stderr,"[1] ./hgm_jack-n \n");
1.28 ! takayama 1651: fprintf(stderr,"[2] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --degree 15 --automatic 0\n");
! 1652: fprintf(stderr,"[3] ./hgm_jack-n --idata Testdata/tmp-idata2.txt --degree 15 >test2.txt\n");
1.14 takayama 1653: fprintf(stderr," ./hgm_w-n --idata test2.txt --gnuplotf test-g\n");
1.1 takayama 1654: fprintf(stderr," gnuplot -persist <test-g-gp.txt\n");
1.28 ! takayama 1655: fprintf(stderr,"[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --automatic 1 --assigned_series_error=1e-12\n");
! 1656: fprintf(stderr,"[5] ./hgm_jack-n --idata Testdata/tmp-idata4.txt\n");
1.13 takayama 1657: return(0);
1.1 takayama 1658: }
1659:
1.10 takayama 1660: static int setParamDefault() {
1.1 takayama 1661: int rank;
1662: int i;
1.3 takayama 1663: Mg = M_n_default ;
1.1 takayama 1664: rank = imypower(2,Mg);
1665: Beta = (double *)mymalloc(sizeof(double)*Mg);
1666: for (i=0; i<Mg; i++) Beta[i] = 1.0+i;
1667: Ng = (double *)mymalloc(sizeof(double)); *Ng = 3.0;
1668: Iv = (double *)mymalloc(sizeof(double)*rank);
1669: Iv2 = (double *)mymalloc(sizeof(double)*rank);
1670: for (i=0; i<rank; i++) Iv[i] = 0;
1671: Ef = 0;
1672: Ef2 = 0.01034957388338225707;
1673: if (M_n == 2) {
1674: Iv2[0] = 1.58693;
1675: Iv2[1] = 0.811369;
1676: Iv2[2] = 0.846874;
1677: Iv2[3] = 0.413438;
1678: }
1679: X0g = (Beta[0]/Beta[Mg-1])*0.5;
1680: Hg = 0.001;
1681: Dp = 1;
1682: Xng = 10.0;
1.13 takayama 1683: return(0);
1.1 takayama 1684: }
1685:
1.10 takayama 1686: static int next(struct SFILE *sfp,char *s,char *msg) {
1.14 takayama 1687: static int check=1;
1688: char *ng="%Ng=";
1689: int i;
1.1 takayama 1690: s[0] = '%';
1691: while (s[0] == '%') {
1.12 takayama 1692: if (!mh_fgets(s,SMAX,sfp)) {
1693: fprintf(stderr,"Data format error at %s\n",msg);
1694: mh_exit(-1);
1695: }
1.14 takayama 1696: if (check && (strncmp(msg,ng,4)==0)) {
1697: if (strncmp(s,ng,4) != 0) {
1698: fprintf(stderr,"Warning, there is no %%Ng= at the border of Beta's and Ng, s=%s\n",s);
1699: }
1700: check=0;
1701: }
1.12 takayama 1702: if (s[0] != '%') return(0);
1.1 takayama 1703: }
1.13 takayama 1704: return(0);
1.1 takayama 1705: }
1.10 takayama 1706: static int setParam(char *fname) {
1.1 takayama 1707: int rank;
1708: char s[SMAX];
1.3 takayama 1709: struct SFILE *fp;
1.1 takayama 1710: int i;
1.19 takayama 1711: struct mh_token tk;
1.1 takayama 1712: if (fname == NULL) return(setParamDefault());
1713:
1714: Sample = 0;
1.3 takayama 1715: if ((fp=mh_fopen(fname,"r",JK_byFile)) == NULL) {
1.12 takayama 1716: if (JK_byFile) fprintf(stderr,"File %s is not found.\n",fp->s);
1717: mh_exit(-1);
1.1 takayama 1718: }
1719: next(fp,s,"Mg(m)");
1720: sscanf(s,"%d",&Mg);
1721: rank = imypower(2,Mg);
1722:
1723: Beta = (double *)mymalloc(sizeof(double)*Mg);
1724: for (i=0; i<Mg; i++) {
1725: next(fp,s,"Beta");
1.12 takayama 1726: sscanf(s,"%lf",&(Beta[i]));
1.1 takayama 1727: }
1728:
1729: Ng = (double *)mymalloc(sizeof(double));
1.14 takayama 1730: next(fp,s,"%Ng= (freedom parameter n)");
1.1 takayama 1731: sscanf(s,"%lf",Ng);
1732:
1733: next(fp,s,"X0g(initial point)");
1734: sscanf(s,"%lf",&X0g);
1.14 takayama 1735:
1.1 takayama 1736: Iv = (double *)mymalloc(sizeof(double)*rank);
1737: for (i=0; i<rank; i++) {
1.12 takayama 1738: next(fp,s,"Iv(initial values)");
1.14 takayama 1739: if (strncmp(s,"*",1)==0) {
1740: for (i=0; i<rank; i++) Iv[i] = 0.0;
1741: break;
1742: }
1.12 takayama 1743: sscanf(s,"%lg",&(Iv[i]));
1.1 takayama 1744: }
1745:
1746: next(fp,s,"Ef(exponential factor)");
1.14 takayama 1747: if (strncmp(s,"*",1)==0) Ef=0.0;
1748: else sscanf(s,"%lg",&Ef);
1.1 takayama 1749:
1750: next(fp,s,"Hg (step size of rk)");
1751: sscanf(s,"%lf",&Hg);
1752:
1753: next(fp,s,"Dp (data sampling period)");
1754: sscanf(s,"%d",&Dp);
1755:
1756: next(fp,s,"Xng (the last point, cf. --xmax)");
1757: sscanf(s,"%lf",&Xng);
1.19 takayama 1758:
1759: /* Reading the optional parameters */
1760: while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) {
1761: /* expect ID */
1762: if (tk.type != MH_TOKEN_ID) {
1763: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1764: }
1765: if (strcmp(s,"automatic")==0) {
1766: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1767: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1768: }
1769: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1770: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1771: }
1772: M_automatic = tk.ival;
1773: continue;
1774: }
1775: if (strcmp(s,"assigned_series_error")==0) {
1776: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1777: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1778: }
1779: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
1780: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1781: }
1782: M_assigned_series_error = tk.dval;
1783: continue;
1784: }
1785: if (strcmp(s,"x0value_min")==0) {
1786: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1787: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1788: }
1789: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
1790: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1791: }
1792: M_x0value_min = tk.dval;
1793: continue;
1794: }
1795: if ((strcmp(s,"Mapprox")==0) || (strcmp(s,"degree")==0)) {
1796: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1797: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1798: }
1799: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1800: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1801: }
1802: Mapprox = tk.ival;
1803: continue;
1804: }
1805: if (strcmp(s,"X0g_bound")==0) {
1806: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1807: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1808: }
1809: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {
1810: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1811: }
1812: M_X0g_bound = tk.dval;
1813: continue;
1814: }
1.25 takayama 1815: if (strcmp(s,"show_autosteps")==0) {
1816: if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
1817: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1818: }
1819: if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
1820: fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
1821: }
1822: M_show_autosteps = tk.ival;
1823: continue;
1824: }
1.19 takayama 1825: fprintf(stderr,"Unknown ID at %s\n",s); mh_exit(-1);
1826: }
1.3 takayama 1827: mh_fclose(fp);
1.13 takayama 1828: return(0);
1.1 takayama 1829: }
1830:
1.10 takayama 1831: static int showParam(struct SFILE *fp,int fd) {
1.1 takayama 1832: int rank,i;
1.3 takayama 1833: char swork[1024];
1.8 takayama 1834: if (fd) {
1835: fp = mh_fopen("stdout","w",1);
1836: }
1.1 takayama 1837: rank = imypower(2,Mg);
1.3 takayama 1838: sprintf(swork,"%%Mg=\n%d\n",Mg); mh_fputs(swork,fp);
1.1 takayama 1839: for (i=0; i<Mg; i++) {
1.12 takayama 1840: sprintf(swork,"%%Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp);
1.1 takayama 1841: }
1.3 takayama 1842: sprintf(swork,"%%Ng=\n%lf\n",*Ng); mh_fputs(swork,fp);
1843: sprintf(swork,"%%X0g=\n%lf\n",X0g); mh_fputs(swork,fp);
1.1 takayama 1844: for (i=0; i<rank; i++) {
1.12 takayama 1845: sprintf(swork,"%%Iv[%d]=\n%lg\n",i,Iv[i]); mh_fputs(swork,fp);
1846: if (Sample && (M_n == 2) && (X0g == 0.3)) {
1847: sprintf(swork,"%%Iv[%d]-Iv2[%d]=%lg\n",i,i,Iv[i]-Iv2[i]); mh_fputs(swork,fp);
1848: }
1.3 takayama 1849: }
1850: sprintf(swork,"%%Ef=\n%lg\n",Ef); mh_fputs(swork,fp);
1851: sprintf(swork,"%%Hg=\n%lf\n",Hg); mh_fputs(swork,fp);
1852: sprintf(swork,"%%Dp=\n%d\n",Dp); mh_fputs(swork,fp);
1853: sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp);
1.19 takayama 1854:
1855: sprintf(swork,"%%%% Optional paramters\n"); mh_fputs(swork,fp);
1856: sprintf(swork,"#success=%d\n",M_mh_t_success); mh_fputs(swork,fp);
1857: sprintf(swork,"#automatic=%d\n",M_automatic); mh_fputs(swork,fp);
1858: sprintf(swork,"#series_error=%lg\n",M_series_error); mh_fputs(swork,fp);
1859: sprintf(swork,"#recommended_abserr\n"); mh_fputs(swork,fp);
1860: sprintf(swork,"%%abserror=%lg\n",M_recommended_abserr); mh_fputs(swork,fp);
1.27 takayama 1861: if (M_recommended_relerr < MH_RELERR_DEFAULT) {
1862: sprintf(swork,"%%relerror=%lg\n",M_recommended_relerr); mh_fputs(swork,fp);
1863: }
1.19 takayama 1864: sprintf(swork,"#mh_t_value=%lg # Value of matrix hg at X0g.\n",M_mh_t_value); mh_fputs(swork,fp);
1865: sprintf(swork,"# M_m=%d # Approximation degree of matrix hg.\n",M_m); mh_fputs(swork,fp);
1866: sprintf(swork,"#beta_i_x_o2_max=%lg #max(|beta[i]*x|/2)\n",M_beta_i_x_o2_max); mh_fputs(swork,fp);
1867: sprintf(swork,"#beta_i_beta_j_min=%lg #min(|beta[i]-beta[j]|)\n",M_beta_i_beta_j_min); mh_fputs(swork,fp);
1.13 takayama 1868: return(0);
1.1 takayama 1869: }
1870:
1.2 takayama 1871: static double gammam(double a,int n) {
1.1 takayama 1872: double v,v2;
1873: int i;
1874: v=mypower(sqrt(M_PI),(n*(n-1))/2);
1875: v2=0;
1876: for (i=1; i<=n; i++) {
1877: v2 += lgamma(a-((double)(i-1))/2.0); /* not for big n */
1878: }
1879: if (Debug) printf("gammam(%lg,%d)=%lg\n",a,n,v*exp(v2));
1880: return(v*exp(v2));
1881: }
1882:
1.2 takayama 1883: static double iv_factor(void) {
1.1 takayama 1884: double v1;
1885: double t;
1886: double b;
1887: double detSigma;
1888: double c;
1889: int i,n;
1890: n = (int) (*Ng);
1891: v1= mypower(sqrt(X0g),n*M_n);
1892: t = 0.0;
1893: for (i=0; i<M_n; i++) t += -X0g*Beta[i];
1894: v1 = v1*exp(t);
1895:
1896: b = 1; for (i=0; i<M_n; i++) b *= Beta[i];
1897: detSigma = 1.0/(b*mypower(2.0,M_n));
1898:
1899: c = gammam(((double)(M_n+1))/2.0, M_n)/
1.12 takayama 1900: ( mypower(sqrt(2), M_n*n)*mypower(sqrt(detSigma),n)*gammam(((double)(M_n+n+1))/2.0,M_n) );
1.1 takayama 1901: return( c*v1);
1902: }
1903:
1904:
1905:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to jack-n.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / hgm / mh / src