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