OpenXM/src/hgm/mh/src/jack-n.c - diff

Return to jack-n.c CVS log

Up to [local] / OpenXM / src / hgm / mh / src

Diff for /OpenXM/src/hgm/mh/src/jack-n.c between version 1.5 and 1.13

version 1.5, 2013/02/21 12:14:08

version 1.13, 2013/03/07 05:23:31

Line 5

#include <string.h>

#include "sfile.h"

$OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.4 2013/02/21 07:51:57 takayama Exp $

$OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.12 2013/03/07 03:00:43 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.

cf. 20-my-note-mh-jack-n.pdf, /Movies/oxvh/priv-j1/2011-12-14-ke-mh-jack.mov

Todo:

1. Cash the transposes of partitions.

2. Use the recurrence to obtain beta().

3. Easier input data file format for mh-n.c

Changelog:

2012.02.21, porting to OpenXM/src/hgm

2011.12.22, --table option, which is experimental.

2011.12.19, bug fix. jjk() should return double. It can become more than max int.

2011.12.15, mh.r --> jack-n.c

/****** from mh-n.c *****/

static int JK_byFile=1;

static int JK_deallocate=0;

#define M_n_default 3

#define Sample_default 1

static int M_n=0;

Line 54 static double Ef2;

Line 55 static double Ef2;

#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 Debug = 1;

static int Debug = 0;

static int Alpha = 2; /* 2 implies the zonal polynomial */

static int *Darray = NULL;

static int **Parray = NULL; /* array of partitions of size M_n */

Line 94 static double M_rel_error=0.0; /* relative errors */

Line 95 static double M_rel_error=0.0; /* relative errors */

/* prototypes */

static void *mymalloc(int size);

static myfree(void *p);

static int myfree(void *p);

static myerror(char *s);

static int myerror(char *s);

static double jack1(int K);

static double jack1diff(int k);

static double xval(int i,int p); /* x_i^p */

Line 103 static int mysum(int L[]);

Line 104 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 test_ptrans();

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[]);

Line 113 static double mypower(double x,int n);

Line 114 static double mypower(double x,int n);

static double qk(int K[],double A[A_LEN],double B[B_LEN]);

static int bb(int N[],int K[],int M[],int I,int J);

static double beta(int K[],int M[]);

static printp(int kappa[]);

static int printp(int kappa[]);

static printp2(int kappa[]);

static int printp2(int kappa[]);

static test_beta();

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();

Line 131 static int pListHS2(int From,int To,int Kap[]);

Line 132 static int pListHS2(int From,int To,int Kap[]);

static void hsExec_0();

static int pmn(int M,int N);

static int *cloneP(int a[]);

static copyP(int p[],int a[]);

static int copyP(int p[],int a[]);

static void pExec_darray(void);

static genDarray2(int M,int N);

static int genDarray2(int M,int N);

static isHStrip(int Kap[],int Nu[]);

static int isHStrip(int Kap[],int Nu[]);

static void hsExec_beta(void);

static genBeta(int Kap[]);

static int genBeta(int Kap[]);

static checkBeta1();

static int checkBeta1();

static int psublen(int Kap[],int Mu[]);

static genJack(int M,int N);

static int genJack(int M,int N);

static checkJack1(int M,int N);

static int checkJack1(int M,int N);

static checkJack2(int M,int N);

static int checkJack2(int M,int N);

static mtest1b();

static int mtest1b();

static int imypower(int x,int n);

static usage();

static int usage();

static setParamDefault();

static int setParamDefault();

static next(struct SFILE *fp,char *s,char *msg);

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);

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);

Line 163 int jk_initializeWorkArea();

Line 164 int jk_initializeWorkArea();

int jk_freeWorkArea() {

/* bug, p in the cloneP will not be deallocated.

Nk in genDarray2 will not be deallocated.

int i;

JK_deallocate=1;

if (Darray) {myfree(Darray); Darray=NULL;}

if (Parray) {myfree(Parray); Parray=NULL;}

if (ParraySize) {myfree(ParraySize); ParraySize=NULL;}

if (M_beta_0) {myfree(M_beta_0); M_beta_0=NULL;}

if (M_beta_1) {myfree(M_beta_1); M_beta_1=NULL;}

if (M_jack) {

for (i=0; M_jack[i] != NULL; i++) {

myfree(M_jack[i]); M_jack[i] = NULL;

if (Debug) printf("Free M_jack[%d]\n",i);

}

myfree(M_jack[i]); M_jack[i] = NULL;

myfree(M_jack); M_jack=NULL;

}

myfree(M_jack); M_jack=NULL;

}

if (M_qk) {myfree(M_qk); M_qk=NULL;}

if (P_pki) {myfree(P_pki); P_pki=NULL;}

JK_deallocate=0;

return(0);

}

int jk_initializeWorkArea() {

int i,j;

JK_deallocate=1;

xval(0,0);

JK_deallocate=0;

Darray=NULL;

Parray=NULL;

ParraySize=NULL;

Line 206 int jk_initializeWorkArea() {

Line 214 int jk_initializeWorkArea() {

Sample = Sample_default;

Xng=0.0;

M_n=0;

return(0);

}

static void *mymalloc(int size) {

void *p;

if (Debug) printf("mymalloc(%d)\n",size);

p = (void *)malloc(size);

p = (void *)mh_malloc(size);

if (p == NULL) {

fprintf(stderr,"No more memory.\n");

mh_exit(-1);

}

return(p);

}

static myfree(void *p) { free(p); }

static int myfree(void *p) {

static myerror(char *s) { fprintf(stderr,"%s: type in control-C\n",s); getchar(); getchar();}

if (Debug) printf("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 double jack1(int K) {

double F;

Line 252 static double xval(int ii,int p) { /* x_i^p */

Line 264 static double xval(int ii,int p) { /* x_i^p */

extern double M_x[];

double F;

int i,j;

static init=0;

static int init=0;

if (JK_deallocate) { init=0; return(0.0);}

if (!init) {

for (i=1; i<=M_n; i++) {

for (j=0; j<M_m_MAX; j++) {

if (j != 0) {

Xarray[i-1][j] = M_x[i-1]*Xarray[i-1][j-1];

}else{

Xarray[i-1][0] = 1;

}

init = 1;

Line 284 static int mysum(int L[]) {

Line 297 static int mysum(int L[]) {

}

(3,2,2,0,0) --> 3

static int plength(int P[]) {

int I;

Line 303 static int plength_t(int P[]) {

Line 316 static int plength_t(int P[]) {

}

ptrans(P) returns Pt

static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m] */

extern int M_m;

Line 318 static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m]

Line 331 static void ptrans(int P[M_nmx],int Pt[]) { /* Pt[M_m]

}

static test_ptrans() {

static int test_ptrans() {

extern int M_m;

int p[M_n0]={5,3,2};

int pt[10];

Line 326 static test_ptrans() {

Line 339 static test_ptrans() {

M_m = 10;

ptrans(p,pt);

if (Debug) {for (i=0; i<10; i++) printf("%d,",pt[i]); printf("\n");}

return(0);

}

Line 418 static double mypower(double x,int n) {

Line 432 static double mypower(double x,int n) {

return(v);

}

/* Q_kappa

static double qk(int K[],double A[A_LEN],double B[B_LEN]) {

extern int Alpha;

int P,Q,I;

Line 438 static double qk(int K[],double A[A_LEN],double B[B_LE

Line 452 static double qk(int K[],double A[A_LEN],double B[B_LE

}

B^nu_{kappa,mu}(i,j)

bb(N,K,M,I,J)

static int bb(int N[],int K[],int M[],int I,int J) {

int Kp[M_m_MAX]; int Mp[M_m_MAX];

Line 447 static int bb(int N[],int K[],int M[],int I,int J) {

Line 461 static int bb(int N[],int K[],int M[],int I,int J) {

ptrans(M,Mp);

printp(K); printf("K<--, "); printp2(Kp); printf("<--Kp\n");

printp(M); printf("M<--, "); printp2(Mp); printf("<--Mp\n");

if ((plength_t(Kp) < J) || (plength_t(Mp) < J)) return(hdk(N,I,J));

Line 484 static double beta(int K[],int M[]) {

Line 498 static double beta(int K[],int M[]) {

return(V);

}

static printp(int kappa[]) {

static int printp(int kappa[]) {

int i;

printf("(");

for (i=0; i<M_n; i++) {

if (i <M_n-1) printf("%d,",kappa[i]);

else printf("%d)",kappa[i]);

}

return(0);

}

static printp2(int kappa[]) {

static int printp2(int kappa[]) {

int i,ell;

printf("(");

ell = plength_t(kappa);

Line 500 static printp2(int kappa[]) {

Line 515 static printp2(int kappa[]) {

if (i <ell+1-1) printf("%d,",kappa[i]);

else printf("%d)",kappa[i]);

}

return(0);

}

static test_beta() {

static int test_beta() {

int kappa[M_n0]={2,1,0};

int mu1[M_n0]={1,0,0};

int mu2[M_n0]={1,1,0};

Line 516 static test_beta() {

Line 532 static test_beta() {

cf. w1m.rr

matrix hypergeometric by jack

N variables, up to degree M.

/* todo

def mhgj(A,B,N,M) {

F = 0;

P = partition_a(N,M);

for (I=0; I<length(P); I++) {

K = P[I];

F += qk(K,A,B)*jack(K,N);

}

return(F);

}

Line 576 static double q3_5(double A[],double B[],int K[],int I

Line 592 static double q3_5(double A[],double B[],int K[],int I

V=1;

for (J=1; J<=P; J++) {

V *= (A[J-1]+C);

}

for (J=1; J<=Q; J++) {

V /= (B[J-1]+C);

}

for (J=1; J<=K[I-1]-1; J++) {

Ej = D-J*Alpha+Kp[J-1];

Gj = Ej+1;

V *= (Gj-Alpha)*Ej/(Gj*(Ej+Alpha));

}

for (J=1; J<=I-1; J++) {

Fj=K[J-1]*Alpha-J-D;

Hj=Fj+Alpha;

Lj=Hj*Fj;

V *= (Lj-Fj)/(Lj+Hj);

}

return(V);

}

Line 621 static void mtest4b() {

Line 637 static void mtest4b() {

/* main() { mtest4(); mtest4b(); } */

/* nk in (4.1),

static int nk(int KK[]) {

extern int *Darray;

int N,I,Ki;

Line 694 static int pListPartition(int M,int N) {

Line 710 static int pListPartition(int M,int N) {

static int pListPartition2(int Less,int From,int To, int M) {

int I,R;

mh_check_intr(100);

if (To < From) {

(*M_pExec)(); return(0);

}

for (I=1; (I<=Less) && (I<=M) ; I++) {

M_kap[From-1] = I;

Line 709 static int pListPartition2(int Less,int From,int To, i

Line 726 static int pListPartition2(int Less,int From,int To, i

static void pExec_0() {

if (Debug) {

printf("M_kap=");

printp(M_kap);

printf("\n");

}

/* Test.

Compare pListPartition(4,3); genDarray(4,3);

Compare pListPartition(5,3); genDarray(5,3);

main() {

M_pExec = pExec_0;

pListPartition(5,3);

}

Line 740 static int pListHS(int Kap[],int N) {

Line 757 static int pListHS(int Kap[],int N) {

HS_n = N;

/* Clear HS_mu. Do not forget when N < M_n */

for (i=0; i<M_n; i++) HS_mu[i] = 0;

pListHS2(1,N,Kap);

return(pListHS2(1,N,Kap));

}

static int pListHS2(int From,int To,int Kap[]) {

Line 763 static void hsExec_0() {

Line 780 static void hsExec_0() {

pListHS([0,4,2,1],3);

main() {

int Kap[3]={4,2,1};

HS_hsExec = hsExec_0;

pListHS(Kap,3);

}

/* The number of partitions <= M, with N parts.

(0,0,...,0) is excluded.

#define aP_pki(i,j) P_pki[(i)*(M+1)+(j)]

static int pmn(int M,int N) {

Line 805 static int pmn(int M,int N) {

Line 822 static int pmn(int M,int N) {

}

main() {pmn(4,3); printf("P_pmn=%d\n",P_pmn);}

static int *cloneP(int a[]) {

Line 815 static int *cloneP(int a[]) {

Line 832 static int *cloneP(int a[]) {

for (i=0; i<M_n; i++) p[i] = a[i];

return(p);

}

static copyP(int p[],int a[]) {

static int copyP(int p[],int a[]) {

int i;

for (i=0; i<M_n; i++) p[i] = a[i];

return(0);

}

static void pExec_darray(void) {

Line 832 static void pExec_darray(void) {

Line 850 static void pExec_darray(void) {

ParraySize[DR_parray] = mysum(K);

DR_parray++;

}

static genDarray2(int M,int N) {

static int genDarray2(int M,int N) {

extern int *Darray;

extern int **Parray;

extern int DR_parray;

Line 860 static genDarray2(int M,int N) {

Line 878 static genDarray2(int M,int N) {

Nk = (int *) mymalloc(sizeof(int)*(Pmn+1));

for (I=0; I<Pmn; I++) Nk[I] = I;

for (I=0; I<Pmn; I++) {

K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */

mh_check_intr(100);

Ksize = plength(K);

K = L[I]; /* N_K = I; D[N_K] = N_(K,1) */

if (Ksize >= M_n) {

Ksize = plength(K);

if (Debug) {fprintf(stderr,"Ksize >= M_n\n");}

if (Ksize >= M_n) {

continue;

if (Debug) {fprintf(stderr,"Ksize >= M_n\n");}

}

continue;

for (i=0; i<M_n; i++) Kone[i] = 0;

}

for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1;

for (i=0; i<M_n; i++) Kone[i] = 0;

for (J=0; J<Pmn; J++) {

for(J=0; J<Ksize; J++) Kone[J]=K[J]; Kone[Ksize] = 1;

if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */

for (J=0; J<Pmn; J++) {

}

if (myeq(L[J],Kone)) Darray[I] = J; /* J is the next of I */

}

if (Debug) {

printf("Darray=\n");

for (i=0; i<Pmn; i++) printf("%d\n",Darray[i]);

printf("-----------\n");

}

return(0);

}

/* main() { genDarray2(4,3);} */

/* M_beta_0[*] value of beta_{kappa,mu}, M_beta_1[*] N_mu */

static isHStrip(int Kap[],int Nu[]) {

static int isHStrip(int Kap[],int Nu[]) {

int N1,N2,I,P;

N1 = plength(Kap); N2 = plength(Nu);

if (N2 > N1) return(0);

Line 925 static void hsExec_beta(void) {

Line 945 static void hsExec_beta(void) {

Done=0;

for (J=M_beta_pt-1; J>=0; J--) {

if (M_beta_1[J] == Nnu) {

K=I+1;

if (Debug) {

printf("Found at J=%d, K=%d, q3_10(Kap,Nu,K)=%lf,Nu,Mu= \n",

J,K,q3_10(M_beta_kap,Nu,K));

printp(Nu); printf("\n");

printp(Mu); printf("\n");

}

/* Check other conditions. See Numata's mail on Dec 24, 2011. */

rrMax = Nu[I]-1;

if ((plength_t(Kapt) < rrMax) || (plength_t(Nut) < rrMax)) {

if (Debug) printf(" is not taken (length). \n");

break;

}

OK=1;

for (RR=0; RR<rrMax; RR++) {

if (Kapt[RR] != Nut[RR]) { OK=0; break;}

}

if (!OK) { if (Debug) printf(" is not taken.\n"); break; }

/* check done. */

M_beta_0[M_beta_pt]=M_beta_0[J]*q3_10(M_beta_kap,Nu,K);

Done = 1; break;

}

if (Done) break; else Nu[I]--;

}

if (!Done) {

if (Debug) printf("BUG: not found M_beta_pt=%d.\n",M_beta_pt);

/* M_beta_0[M_beta_pt] = NAN; /* error("Not found."); */

/* M_beta_0[M_beta_pt] = NAN; error("Not found."); */

M_beta_0[M_beta_pt] = beta(M_beta_kap,Mu);

}

/* Fix the bug of mh.rr */

M_beta_pt++;

}

static genBeta(int Kap[]) {

static int genBeta(int Kap[]) {

extern double *M_beta_0;

extern int *M_beta_1;

extern int M_beta_pt;

Line 968 static genBeta(int Kap[]) {

Line 988 static genBeta(int Kap[]) {

if (Debug) {printp(Kap); printf("<-Kappa, P_pmn=%d\n",P_pmn);}

/* M_beta = newmat(2,P_pmn+1); */

M_beta_0 = (double *)mymalloc(sizeof(double)*(P_pmn+1));

M_beta_1 = (int *)mymalloc(sizeof(double)*(P_pmn+1));

M_beta_1 = (int *)mymalloc(sizeof(int)*(P_pmn+1));

M_beta_pt = 0;

for (I=0; I<=P_pmn; I++) {M_beta_0[I] = NAN; M_beta_1[I] = -1;}

N = plength(Kap);

Line 982 static genBeta(int Kap[]) {

Line 1002 static genBeta(int Kap[]) {

genBeta([2,1,1]);

static checkBeta1() {

static int checkBeta1() {

int Kap[3] = {2,2,0};

int Kap2[3] = {2,1,0};

int I;

Line 993 static checkBeta1() {

Line 1013 static checkBeta1() {

for (I=0; I<M_beta_pt; I++) {

Mu = Parray[M_beta_1[I]];

Beta_km = M_beta_0[I];

if (Debug) {

printp(Kap); printf("<--Kap, ");

printp(Mu); printf("<--Mu,");

printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap,Mu));

}

if (Debug) printf("-------------------------------------\n");

genBeta(Kap2);

for (I=0; I<M_beta_pt; I++) {

Mu = Parray[M_beta_1[I]];

Beta_km = M_beta_0[I];

if (Debug) {

printp(Kap2); printf("<--Kap, ");

printp(Mu); printf("<--Mu,");

printf("Beta_km(by table)=%lf, beta(Kap,Mu)=%lf\n",Beta_km,beta(Kap2,Mu));

}

return(0);

}

def checkBeta2() {

genDarray2(3,3);

Kap = [2,1,0];

printf("Kap=%a\n",Kap);

genBeta(Kap);

for (I=0; I<M_beta_pt; I++) {

Mu = Parray[M_beta[1][I]];

Beta_km = M_beta[0][I];

printf("Mu=%a,",Mu);

printf("Beta_km(by table)=%a, beta(Kap,Mu)=%a\n",Beta_km,beta(Kap,Mu));

}

/* main() { checkBeta1(); } */

Line 1045 static int psublen(int Kap[],int Mu[]) {

Line 1066 static int psublen(int Kap[],int Mu[]) {

/* Table of Jack polynomials

Jack[1]* one variable.

Jack[2]* two variables.

...

Jack[M_n]* n variables.

Jack[P][J]*

D^J(P variables jack of p variables). Example. J=001 d_1, 010 d_2, 100 d_3

0<=J<=2^{M_n}-1

Jack[P][J][nk(Kappa)] Jack_Kappa, Kappa is a partition.

0<=nk(Kappa)<=pmn(M_m,M_n)

#define aM_jack(i,j,k) ((M_jack[i])[(j)*(Pmn+1)+(k)])

static genJack(int M,int N) {

static int genJack(int M,int N) {

extern double **M_jack;

extern int M_2n;

extern int P_pmn;

Line 1095 static genJack(int M,int N) {

Line 1116 static genJack(int M,int N) {

if (M_df) {

for (J=1; J<M_2n; J++) {

if (J >= 2^I) aM_jack(I,J,K) = 0;

else aM_jack(I,J,K) = NAN;

}

Line 1120 static genJack(int M,int N) {

Line 1141 static genJack(int M,int N) {

for (H=0; H<M_beta_pt; H++) {

Nk = M_beta_1[H];

Mu = Parray[Nk];

if (UseTable) {

Beta_km = M_beta_0[H];

}else{

Beta_km = beta(Kap,Mu);

/* do not use the M_beta table. It's buggy. UseTable is experimental.*/

}

if (Debug) {printf("Nv(number of variables)=%d, Beta_km=%lf, Mu=",Nv,Beta_km);

printp(Mu); printf("\n");}

P = psublen(Kap,Mu);

Jack += aM_jack(Nv-1,0,Nk)*Beta_km*xval(Nv,P); /* util_v(x,[Nv])^P;*/

if (Debug) printf("xval(%d,%d)=%lf\n",Nv,P,xval(Nv,P));

}

aM_jack(Nv,0,K) = Jack;

if (M_df) {

/* The case of M_df > 0. */

for (J=1; J<M_2n; J++) {

Jack = 0;

mh_check_intr(100);

for (H=0; H<M_beta_pt; H++) {

Jack = 0;

Nk = M_beta_1[H];

for (H=0; H<M_beta_pt; H++) {

Mu = Parray[Nk];

Nk = M_beta_1[H];

if (UseTable) {

Mu = Parray[Nk];

Beta_km = M_beta_0[H];

if (UseTable) {

}else{

Beta_km = M_beta_0[H];

Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */

}else{

}

Beta_km = beta(Kap,Mu); /* do not use the M_beta table. It's buggy. */

if (Debug) {printf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);

}

printp(Mu); printf("\n"); }

if (Debug) {printf("M_df: Nv(number of variables)=%d, Beta_km=%lf, Mu= ",Nv,Beta_km);

P = psublen(Kap,Mu);

printp(Mu); printf("\n"); }

if (J & (1 << (Nv-1))) {

P = psublen(Kap,Mu);

JJ = J & ((1 << (Nv-1)) ^ 0xffff); /* NOTE!! Up to 16 bits. mh-15 */

if (J & (1 << (Nv-1))) {

if (P != 0) {

JJ = J & ((1 << (Nv-1)) ^ 0xffff); /* NOTE!! Up to 16 bits. mh-15 */

Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1);

if (P != 0) {

}

Jack += aM_jack(Nv-1,JJ,Nk)*Beta_km*P*xval(Nv,P-1);

}else{

}

Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P);

}else{

}

Jack += aM_jack(Nv-1,J,Nk)*Beta_km*xval(Nv,P);

}

aM_jack(Nv,J,K) = Jack;

}

if (Debug) printf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);

aM_jack(Nv,J,K) = Jack;

if (Debug) printf("aM_jack(%d,%d,%d) = %lf\n",Nv,J,K,Jack);

} /* end of J loop */

}

Line 1166 static genJack(int M,int N) {

Line 1189 static genJack(int M,int N) {

/* checkJack1(3,3)

static checkJack1(int M,int N) {

static int checkJack1(int M,int N) {

int I,K;

extern int P_pmn;

extern double M_x[];

Line 1186 static checkJack1(int M,int N) {

Line 1209 static checkJack1(int M,int N) {

}

for (I=1; I<=N; I++) printf("%lf, ",M_x[I-1]);

printf("<--x\n");

return(0);

}

/*main() { checkJack1(3,3); }*/

static checkJack2(int M,int N) {

static int checkJack2(int M,int N) {

int I,K,J;

extern int P_pmn;

extern double M_x[];

extern M_df;

extern int M_df;

int Pmn; /* used in aM_jack */

M_df=1;

/* initialize x vars. */

Line 1215 static checkJack2(int M,int N) {

Line 1239 static checkJack2(int M,int N) {

for (I=1; I<=N; I++) {

for (K=0; K<=P_pmn; K++) {

for (J=0; J<M_2n; J++) {

printp(Parray[K]);

printf("<--Kap, Nv=%d,J(diff)=%d, D^J Jack=%lf\n",

I,J,aM_jack(I,J,K));

}

return(0);

}

/* main() { checkJack2(3,3); } */

Line 1238 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

Line 1263 double mh_t(double A[A_LEN],double B[B_LEN],int N,int

F = 0; F2=0;

M_df=1;

genJack(M,N);

M_qk = (double *)mymalloc(sizeof(double)*P_pmn);

M_qk = (double *)mymalloc(sizeof(double)*(P_pmn+1)); /* found a bug. */

Pmn = P_pmn;

size = ParraySize[P_pmn];

for (K=0; K<=P_pmn; K++) {

Kap = Parray[K];

mh_check_intr(100);

M_qk[K] = qk(Kap,A,B);

Kap = Parray[K];

F += M_qk[K]*aM_jack(N,0,K);

M_qk[K] = qk(Kap,A,B);

if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K);

F += M_qk[K]*aM_jack(N,0,K);

if (Debug) printf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size);

if (ParraySize[K] < size) F2 += M_qk[K]*aM_jack(N,0,K);

if (Debug && (ParraySize[K] == size)) printf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K));

if (Debug) printf("ParraySize[K] = %d, size=%d\n",ParraySize[K],size);

if (Debug && (ParraySize[K] == size)) printf("M_qk[K]=%lg, aM_jack=%lg\n",M_qk[K],aM_jack(N,0,K));

}

M_rel_error = F-F2;

return(F);

Line 1262 double mh_t2(int J) {

Line 1288 double mh_t2(int J) {

F = 0;

Pmn = P_pmn;

for (K=0; K<P_pmn; K++) {

F += M_qk[K]*aM_jack(M_n,J,K);

}

return(F);

}

static mtest1b() {

static int mtest1b() {

double A[1] = {1.5};

double B[1] = {1.5+5};

int I,N,M,J;

Line 1278 static mtest1b() {

Line 1304 static mtest1b() {

}

mh_t(A,B,N,M);

for (J=0; J<M_2n; J++) {

F=mh_t2(J);

printf("J=%d, D^J mh_t=%lf\n",J,F);

}

Line 1307 static int imypower(int x,int n) {

Line 1333 static int imypower(int x,int n) {

return(v);

}

#ifdef STANDALONE

#ifdef STANDALONE2

main(int argc,char *argv[]) {

mh_exit(MH_RESET_EXIT);

/* jk_main(argc,argv);

printf("second run.\n"); */

jk_main(argc,argv);

}

#endif

Line 1334 struct MH_RESULT *jk_main(int argc,char *argv[]) {

Line 1360 struct MH_RESULT *jk_main(int argc,char *argv[]) {

UseTable = 1;

Mapprox=6;

for (i=1; i<argc; i++) {

if (strcmp(argv[i],"--idata")==0) {

inci(i);

setParam(argv[i]); idata=1;

}else if (strcmp(argv[i],"--degree")==0) {

inci(i);

sscanf(argv[i],"%d",&Mapprox);

}else if (strcmp(argv[i],"--x0")==0) {

inci(i);

sscanf(argv[i],"%lg",&X0g);

}else if (strcmp(argv[i],"--notable")==0) {

UseTable = 0;

}else if (strcmp(argv[i],"--debug")==0) {

Debug = 1;

}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);}

JK_byFile = 0;

}else {

fprintf(stderr,"Unknown option %s\n",argv[i]);

usage();

return(NULL);

}

if (!idata) setParam(NULL);

Line 1371 struct MH_RESULT *jk_main(int argc,char *argv[]) {

Line 1397 struct MH_RESULT *jk_main(int argc,char *argv[]) {

sprintf(swork,"%%%%Use --help option to see the help.\n"); mh_fputs(swork,ofp);

sprintf(swork,"%%%%Mapprox=%d\n",Mapprox); mh_fputs(swork,ofp);

if (M_n != Mg) {

myerror("Mg must be equal to M_n\n"); mh_exit(-1);

}

if (Debug) showParam(ofp);

if (Debug) showParam(NULL,1);

for (i=0; i<M_n; i++) {

M_x[i] = Beta[i]*X0g;

}

a[0] = ((double)Mg+1.0)/2.0;

b[0] = ((double)Mg+1.0)/2.0 + ((double) (*Ng))/2.0; /* bug, double check */

if (Debug) printf("Calling mh_t with ([%lf],[%lf],%d,%d)\n",a[0],b[0],M_n,Mapprox);

mh_t(a,b,M_n,Mapprox);

if (imypower(2,M_n) != M_2n) {

sprintf(swork,"M_n=%d,M_2n=%d\n",M_n,M_2n); mh_fputs(swork,ofp);

myerror("2^M_n != M_2n\n"); mh_exit(-1);

}

sprintf(swork,"%%%%M_rel_error=%lg\n",M_rel_error); mh_fputs(swork,ofp);

for (j=0; j<M_2n; j++) {

Iv[j] = mh_t2(j);

}

Ef = iv_factor();

showParam(ofp);

showParam(ofp,0);

/* return the result */

if (!JK_byFile) {

ans->x = X0g;

ans->rank = imypower(2,Mg);

ans->y = (double *)mymalloc(sizeof(double)*(ans->rank));

for (i=0; i<ans->rank; i++) (ans->y)[i] = Iv[i];

ans->size=1;

ans->sfpp = (struct SFILE **)mymalloc(sizeof(struct SFILE *)*(ans->size));

(ans->sfpp)[0] = ofp;

}

if (Debug) printf("jk_freeWorkArea() starts\n");

jk_freeWorkArea();

if (Debug) printf("jk_freeWorkArea() has finished.\n");

return(ans);

}

static usage() {

static int usage() {

fprintf(stderr,"Usages:\n");

fprintf(stderr,"mh-m [--idata input_data_file --x0 x0 --degree approxm]\n");

fprintf(stderr,"\nThe command mh-m [options] generates an input for w-m, Pr({y | y<xmax}), which is the cumulative distribution function of the largest root of the m by m Wishart matrix with n degrees of freedom and the covariantce matrix sigma.\n");

Line 1431 static usage() {

Line 1459 static usage() {

fprintf(stderr,"[6] ./mh-3 --degree 15 >test2.txt\n");

fprintf(stderr," ./w-3 --idata test2.txt --gnuplotf test-g\n");

fprintf(stderr," gnuplot -persist <test-g-gp.txt\n");

return(0);

}

static setParamDefault() {

static int setParamDefault() {

int rank;

int i;

Mg = M_n_default ;

Line 1456 static setParamDefault() {

Line 1485 static setParamDefault() {

Hg = 0.001;

Dp = 1;

Xng = 10.0;

return(0);

}

static next(struct SFILE *sfp,char *s,char *msg) {

static int next(struct SFILE *sfp,char *s,char *msg) {

s[0] = '%';

while (s[0] == '%') {

if (!mh_fgets(s,SMAX,sfp)) {

fprintf(stderr,"Data format error at %s\n",msg);

mh_exit(-1);

}

if (s[0] != '%') return(0);

}

return(0);

}

static setParam(char *fname) {

static int setParam(char *fname) {

int rank;

char s[SMAX];

struct SFILE *fp;

Line 1477 static setParam(char *fname) {

Line 1508 static setParam(char *fname) {

Sample = 0;

if ((fp=mh_fopen(fname,"r",JK_byFile)) == NULL) {

if (JK_byFile) fprintf(stderr,"File %s is not found.\n",fp->s);

mh_exit(-1);

}

next(fp,s,"Mg(m)");

sscanf(s,"%d",&Mg);

Line 1487 static setParam(char *fname) {

Line 1518 static setParam(char *fname) {

Beta = (double *)mymalloc(sizeof(double)*Mg);

for (i=0; i<Mg; i++) {

next(fp,s,"Beta");

sscanf(s,"%lf",&(Beta[i]));

}

Ng = (double *)mymalloc(sizeof(double));

Line 1499 static setParam(char *fname) {

Line 1530 static setParam(char *fname) {

Iv = (double *)mymalloc(sizeof(double)*rank);

for (i=0; i<rank; i++) {

next(fp,s,"Iv(initial values)");

sscanf(s,"%lg",&(Iv[i]));

}

next(fp,s,"Ef(exponential factor)");

Line 1515 static setParam(char *fname) {

Line 1546 static setParam(char *fname) {

next(fp,s,"Xng (the last point, cf. --xmax)");

sscanf(s,"%lf",&Xng);

mh_fclose(fp);

return(0);

}

static showParam(struct SFILE *fp) {

static int showParam(struct SFILE *fp,int fd) {

int rank,i;

char swork[1024];

if (fd) {

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<Mg; i++) {

sprintf(swork,"%%Beta[%d]=\n%lf\n",i,Beta[i]); mh_fputs(swork,fp);

}

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<rank; i++) {

sprintf(swork,"%%Iv[%d]=\n%lg\n",i,Iv[i]); mh_fputs(swork,fp);

if (Sample && (M_n == 2) && (X0g == 0.3)) {

sprintf(swork,"%%Iv[%d]-Iv2[%d]=%lg\n",i,i,Iv[i]-Iv2[i]); mh_fputs(swork,fp);

}

sprintf(swork,"%%Ef=\n%lg\n",Ef); mh_fputs(swork,fp);

sprintf(swork,"%%Hg=\n%lf\n",Hg); mh_fputs(swork,fp);

sprintf(swork,"%%Dp=\n%d\n",Dp); mh_fputs(swork,fp);

sprintf(swork,"%%Xng=\n%lf\n",Xng);mh_fputs(swork,fp);

return(0);

}

static double gammam(double a,int n) {

Line 1568 static double iv_factor(void) {

Line 1604 static double iv_factor(void) {

detSigma = 1.0/(b*mypower(2.0,M_n));

c = gammam(((double)(M_n+1))/2.0, M_n)/

( mypower(sqrt(2), M_n*n)*mypower(sqrt(detSigma),n)*gammam(((double)(M_n+n+1))/2.0,M_n) );

return( c*v1);

}

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>