Annotation of OpenXM/src/hgm/mh/src/mh.c, Revision 1.13
1.13 ! takayama 1: /* $OpenXM: OpenXM/src/hgm/mh/src/mh.c,v 1.12 2013/03/08 04:54:01 takayama Exp $ */
1.1 takayama 2: #include <stdio.h>
3: #include "sfile.h"
4: #include "mh.h"
5: #define WSIZE 1024
1.4 takayama 6: extern int MH_DEBUG;
1.13 ! takayama 7: extern int M_show_autosteps;
1.10 takayama 8: static int imypower(int x,int n) {
1.1 takayama 9: int a,i;
10: a = 1;
11: for (i=0; i<n; i++) a = a*x;
12: return(a);
13: }
14:
15: struct cWishart *new_cWishart(int rank) {
16: struct cWishart *cwp;
17: cwp = (struct cWishart *)mh_malloc(sizeof(struct cWishart));
18: cwp->x=0;
19: cwp->rank=rank;
20: cwp->f = (double *)mh_malloc(sizeof(double)*rank);
21: cwp->aux = cwp->aux2 = NULL;
22: return(cwp);
23: }
24:
1.3 takayama 25: struct cWishart *mh_cwishart_gen(int m,int n,double beta[],double x0,
1.13 ! takayama 26: int approxDeg,double h, int dp, double x,int modep[],
! 27: int automatic,double assigned_series_error,
! 28: int verbose)
! 29: {
1.6 takayama 30: /*
31: modep[0]. Do Koev-Edelman (ignored for now).
32: modep[1]. Do the HGM
33: modep[2]. Return modep[2] intermediate data.
34: */
1.1 takayama 35: struct SFILE *fp;
36: char swork[WSIZE];
37: char *argv[WSIZE];
1.13 ! takayama 38: int argc;
1.1 takayama 39: int i,rank;
40: char *comm;
41: struct MH_RESULT *rp;
42: struct SFILE *sfp;
43: struct cWishart *cw;
44: argv[0]="dummy";
45: argv[1] = "--bystring";
46: argv[2] = "--idata";
47: fp = mh_fopen("","w",0);
48: mh_fputs("%%Mg=\n",fp);
49: sprintf(swork,"%d\n",m); mh_fputs(swork,fp);
50: rank = imypower(2,m);
51: mh_fputs("%%Beta\n",fp);
52: for (i=0; i<m; i++) {
53: sprintf(swork,"%lf\n",beta[i]); mh_fputs(swork,fp);
54: }
55: mh_fputs("%%Ng=\n",fp); /* freedom param */
56: sprintf(swork,"%d\n",n); mh_fputs(swork,fp);
57: mh_fputs("%%X0g=\n",fp); /* initial point */
58: sprintf(swork,"%lf\n",x0); mh_fputs(swork,fp);
59: mh_fputs("%%Iv\n",fp); /* initial values, dummy */
60: for (i=0; i<rank; i++) mh_fputs("0.0\n",fp);
61: mh_fputs("%%Ef=\n1.0\n",fp); /* Below are dummy values */
62: if (h <= 0.0) {fprintf(stderr,"h<=0.0, set to 0.1\n"); h=0.1;}
63: mh_fputs("%%Hg=\n",fp);
64: sprintf(swork,"%lf\n",h); mh_fputs(swork,fp);
65: if (dp < 1) {fprintf(stderr,"dp<1, set to 1\n"); dp=1;}
66: mh_fputs("%%Dp=\n",fp);
67: sprintf(swork,"%d\n",dp); mh_fputs(swork,fp);
68: if (x <= x0) {fprintf(stderr,"x <= x0, set to x=x0+10\n"); x=x0+10;}
69: mh_fputs("%%Xng=\n",fp);
70: sprintf(swork,"%lf\n",x); mh_fputs(swork,fp);
71:
1.13 ! takayama 72: sprintf(swork,"%%automatic=%d\n",automatic); mh_fputs(swork,fp);
! 73: sprintf(swork,"%%assigned_series_error=%lg\n",assigned_series_error); mh_fputs(swork,fp);
! 74: sprintf(swork,"%%show_autosteps=%d\n",verbose); mh_fputs(swork,fp);
! 75:
1.1 takayama 76: comm = (char *)mh_malloc(fp->len +1);
77: mh_outstr(comm,fp->len+1,fp);
78: mh_fclose(fp);
79: argv[3] = comm;
80:
81: argv[4] = "--degree";
82: argv[5] = (char *)mh_malloc(128);
83: sprintf(argv[5],"%d",approxDeg);
84:
1.4 takayama 85: rp=jk_main(6,argv);
1.1 takayama 86: if (rp == NULL) {
87: fprintf(stderr,"rp is NULL.\n"); return(NULL);
88: }
89: cw = new_cWishart(rank);
90: cw->x = rp->x;
91: cw->rank = rp->rank;
1.5 takayama 92: if (rank != cw->rank) {
93: fprintf(stderr,"Rank error.\n"); return(NULL);
94: }
1.4 takayama 95: for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i];
96: sfp = (rp->sfpp)[0];
1.5 takayama 97: cw->aux = (char *) mh_malloc((sfp->len)+1);
98: mh_outstr((char *)(cw->aux),(sfp->len)+1,sfp);
99: /* todo, the following line seems to cause seg fault. */
1.1 takayama 100: /* deallocate the memory */
101: for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]);
102: /* todo, mh_free_??(rp); free Iv's */
1.6 takayama 103: if (!modep[1]) return(cw);
1.1 takayama 104:
1.5 takayama 105: if (MH_DEBUG) printf("\n\n%s\n",(char *)cw->aux);
106: /* This output is strange. */
1.1 takayama 107: /* Starting HGM */
108: argv[3] = (char *)cw->aux;
109: argv[4] = "--dataf";
110: argv[5] = "dummy-dataf";
1.13 ! takayama 111: argc=6;
! 112: if (verbose) {
! 113: argv[6] = "--verbose"; ++argc;
! 114: }
! 115: rp = mh_main(argc,argv);
1.1 takayama 116: if (rp == NULL) {
117: fprintf(stderr,"rp is NULL in the second step.\n"); return(NULL);
118: }
119: cw = new_cWishart(rank);
120: cw->x = rp->x;
121: cw->rank = rp->rank;
122: for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i];
123: sfp = (rp->sfpp)[0];
124: if (sfp) {
125: cw->aux = (char *) mh_malloc(sfp->len+1);
126: mh_outstr((char *)cw->aux,sfp->len+1,sfp);
127: }
128: sfp = (rp->sfpp)[1];
129: if (sfp) {
130: cw->aux2 = (char *) mh_malloc(sfp->len+1);
131: mh_outstr((char *)cw->aux2,sfp->len+1,sfp);
132: }
133: /* deallocate the memory */
134: for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]);
135: mh_freeWorkArea();
136: return(cw);
137: }
138: /* Cumulative probability distribution function of the first eigenvalue of
139: Wishart matrix by Series */
140: struct cWishart *mh_cwishart_s(int m,int n,double beta[],double x0,
1.13 ! takayama 141: int approxDeg,double h, int dp, double x,
! 142: int automatic,double assigned_series_error,
! 143: int verbose)
! 144: {
1.6 takayama 145: int modep[]={1,0,0};
1.13 ! takayama 146: return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep,automatic,assigned_series_error,verbose));
1.1 takayama 147: }
148:
149: /* Cumulative probability distribution function of the first eigenvalue of
150: Wishart matrix by HGM */
151: struct cWishart *mh_cwishart_hgm(int m,int n,double beta[],double x0,
1.13 ! takayama 152: int approxDeg, double h, int dp , double x,
! 153: int automatic,double assigned_series_error,
! 154: int verbose)
1.1 takayama 155: {
1.6 takayama 156: int modep[]={1,1,0};
1.13 ! takayama 157: return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep,automatic,assigned_series_error,verbose));
1.1 takayama 158: }
159:
160: #ifdef STANDALONE
1.12 takayama 161: main(int argc,char *argv[]) {
1.1 takayama 162: double beta[5]={1.0,2.0,3.0,4.0,5.0};
163: struct cWishart *cw;
1.7 takayama 164: struct SFILE *sfp;
165: char *s;
166: char str[1024];
167: double x;
1.12 takayama 168: int i,show;
1.13 ! takayama 169: int strategy=1;
1.12 takayama 170: double err[2]={-1.0,-1.0};
1.13 ! takayama 171: int verbose=0;
1.12 takayama 172: show=1;
173: for (i=1; i<argc; i++) {
174: if (strcmp(argv[i],"--strategy")==0) {
175: i++; sscanf(argv[i],"%d",&strategy);
176: }else if (strcmp(argv[i],"--abserr")==0) {
177: i++; sscanf(argv[i],"%lg",&(err[0]));
178: }else if (strcmp(argv[i],"--relerr")==0) {
179: i++; sscanf(argv[i],"%lg",&(err[1]));
180: }else if (strcmp(argv[i],"--quiet")==0) {
181: show=0;
1.13 ! takayama 182: }else if (strcmp(argv[i],"--verbose")==0) {
! 183: verbose=1;
1.12 takayama 184: }else{
185: fprintf(stderr,"Unknown option.\n");
186: }
187: }
188: mh_set_strategy(strategy,err);
1.13 ! takayama 189: cw=mh_cwishart_hgm(3,5,beta,0.3,7, 0.01,1,10, 1,1e-7,verbose);
1.1 takayama 190: if (cw != NULL) {
191: printf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
192: /* printf("%s",(char *)cw->aux); */
1.11 takayama 193: }
1.13 ! takayama 194: cw=mh_cwishart_hgm(4,5,beta,0.3,7, 0.01,1,10, 1,1e-7,verbose);
1.1 takayama 195: if (cw != NULL) {
196: printf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
1.7 takayama 197: s = (char *)cw->aux;
1.1 takayama 198: /* printf("%s",(char *)cw->aux); */
1.12 takayama 199: if (show) {
200: sfp = mh_fopen(s,"r",0);
201: while (mh_fgets(str,1024,sfp)) {
202: sscanf(str,"%lg",&x); printf("%lg\n",x);
203: }
204: mh_fclose(sfp);
1.7 takayama 205: }
1.1 takayama 206: }
207: }
208: main1() {
209: double beta[5]={1.0,2.0,3.0,4.0,5.0};
210: struct cWishart *cw;
1.13 ! takayama 211: int verbose=1;
! 212: cw=mh_cwishart_s(3,5,beta,0.3,7, 0,0,0, 1,1e-7,verbose);
1.1 takayama 213: if (cw != NULL) {
214: printf("%s",(char *)cw->aux);
215: }
1.13 ! takayama 216: cw=mh_cwishart_s(4,5,beta,0.3,7, 0,0,0, 1,1e-7,verbose);
1.1 takayama 217: if (cw != NULL) {
218: printf("%s",(char *)cw->aux);
219: }
1.13 ! takayama 220: cw=mh_cwishart_s(5,5,beta,0.3,7, 0,0,0, 1,1e-7,verbose);
1.1 takayama 221: if (cw != NULL) {
222: printf("%s",(char *)cw->aux);
223: }
224: }
225: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mh.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