Annotation of OpenXM/src/hgm/mh/src/mh.c, Revision 1.19
1.19 ! takayama 1: /* $OpenXM: OpenXM/src/hgm/mh/src/mh.c,v 1.18 2016/02/15 07:42:07 takayama Exp $ */
1.1 takayama 2: #include <stdio.h>
1.15 takayama 3: #include <string.h>
1.1 takayama 4: #include "sfile.h"
5: #include "mh.h"
6: #define WSIZE 1024
1.16 takayama 7: int MH_DEBUG2=0;
1.4 takayama 8: extern int MH_DEBUG;
1.13 takayama 9: extern int M_show_autosteps;
1.10 takayama 10: static int imypower(int x,int n) {
1.1 takayama 11: int a,i;
12: a = 1;
13: for (i=0; i<n; i++) a = a*x;
14: return(a);
15: }
16:
17: struct cWishart *new_cWishart(int rank) {
18: struct cWishart *cwp;
19: cwp = (struct cWishart *)mh_malloc(sizeof(struct cWishart));
20: cwp->x=0;
21: cwp->rank=rank;
22: cwp->f = (double *)mh_malloc(sizeof(double)*rank);
23: cwp->aux = cwp->aux2 = NULL;
24: return(cwp);
25: }
26:
1.3 takayama 27: struct cWishart *mh_cwishart_gen(int m,int n,double beta[],double x0,
1.13 takayama 28: int approxDeg,double h, int dp, double x,int modep[],
29: int automatic,double assigned_series_error,
30: int verbose)
31: {
1.6 takayama 32: /*
33: modep[0]. Do Koev-Edelman (ignored for now).
34: modep[1]. Do the HGM
35: modep[2]. Return modep[2] intermediate data.
36: */
1.1 takayama 37: struct SFILE *fp;
38: char swork[WSIZE];
39: char *argv[WSIZE];
1.13 takayama 40: int argc;
1.1 takayama 41: int i,rank;
42: char *comm;
43: struct MH_RESULT *rp;
44: struct SFILE *sfp;
45: struct cWishart *cw;
46: argv[0]="dummy";
47: argv[1] = "--bystring";
48: argv[2] = "--idata";
49: fp = mh_fopen("","w",0);
50: mh_fputs("%%Mg=\n",fp);
51: sprintf(swork,"%d\n",m); mh_fputs(swork,fp);
52: rank = imypower(2,m);
53: mh_fputs("%%Beta\n",fp);
54: for (i=0; i<m; i++) {
55: sprintf(swork,"%lf\n",beta[i]); mh_fputs(swork,fp);
56: }
57: mh_fputs("%%Ng=\n",fp); /* freedom param */
58: sprintf(swork,"%d\n",n); mh_fputs(swork,fp);
59: mh_fputs("%%X0g=\n",fp); /* initial point */
60: sprintf(swork,"%lf\n",x0); mh_fputs(swork,fp);
61: mh_fputs("%%Iv\n",fp); /* initial values, dummy */
62: for (i=0; i<rank; i++) mh_fputs("0.0\n",fp);
63: mh_fputs("%%Ef=\n1.0\n",fp); /* Below are dummy values */
1.14 takayama 64: if (h <= 0.0) {oxprintfe("h<=0.0, set to 0.1\n"); h=0.1;}
1.1 takayama 65: mh_fputs("%%Hg=\n",fp);
66: sprintf(swork,"%lf\n",h); mh_fputs(swork,fp);
67: mh_fputs("%%Dp=\n",fp);
68: sprintf(swork,"%d\n",dp); mh_fputs(swork,fp);
1.14 takayama 69: if (x <= x0) {oxprintfe("x <= x0, set to x=x0+10\n"); x=x0+10;}
1.1 takayama 70: mh_fputs("%%Xng=\n",fp);
71: sprintf(swork,"%lf\n",x); mh_fputs(swork,fp);
72:
1.13 takayama 73: sprintf(swork,"%%automatic=%d\n",automatic); mh_fputs(swork,fp);
74: sprintf(swork,"%%assigned_series_error=%lg\n",assigned_series_error); mh_fputs(swork,fp);
75: sprintf(swork,"%%show_autosteps=%d\n",verbose); mh_fputs(swork,fp);
76:
1.1 takayama 77: comm = (char *)mh_malloc(fp->len +1);
78: mh_outstr(comm,fp->len+1,fp);
79: mh_fclose(fp);
80: argv[3] = comm;
81:
82: argv[4] = "--degree";
83: argv[5] = (char *)mh_malloc(128);
84: sprintf(argv[5],"%d",approxDeg);
85:
1.4 takayama 86: rp=jk_main(6,argv);
1.1 takayama 87: if (rp == NULL) {
1.14 takayama 88: oxprintfe("rp is NULL.\n"); return(NULL);
1.1 takayama 89: }
90: cw = new_cWishart(rank);
91: cw->x = rp->x;
92: cw->rank = rp->rank;
1.5 takayama 93: if (rank != cw->rank) {
1.14 takayama 94: oxprintfe("Rank error.\n"); return(NULL);
1.5 takayama 95: }
1.4 takayama 96: for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i];
97: sfp = (rp->sfpp)[0];
1.5 takayama 98: cw->aux = (char *) mh_malloc((sfp->len)+1);
99: mh_outstr((char *)(cw->aux),(sfp->len)+1,sfp);
100: /* todo, the following line seems to cause seg fault. */
1.1 takayama 101: /* deallocate the memory */
1.18 takayama 102: //for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]);
1.1 takayama 103: /* todo, mh_free_??(rp); free Iv's */
1.6 takayama 104: if (!modep[1]) return(cw);
1.1 takayama 105:
1.14 takayama 106: if (MH_DEBUG) oxprintf("\n\n%s\n",(char *)cw->aux);
1.5 takayama 107: /* This output is strange. */
1.1 takayama 108: /* Starting HGM */
109: argv[3] = (char *)cw->aux;
110: argv[4] = "--dataf";
111: argv[5] = "dummy-dataf";
1.13 takayama 112: argc=6;
113: if (verbose) {
114: argv[6] = "--verbose"; ++argc;
115: }
116: rp = mh_main(argc,argv);
1.1 takayama 117: if (rp == NULL) {
1.14 takayama 118: oxprintfe("rp is NULL in the second step.\n"); return(NULL);
1.1 takayama 119: }
120: cw = new_cWishart(rank);
121: cw->x = rp->x;
122: cw->rank = rp->rank;
123: for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i];
124: sfp = (rp->sfpp)[0];
125: if (sfp) {
126: cw->aux = (char *) mh_malloc(sfp->len+1);
127: mh_outstr((char *)cw->aux,sfp->len+1,sfp);
128: }
129: sfp = (rp->sfpp)[1];
130: if (sfp) {
131: cw->aux2 = (char *) mh_malloc(sfp->len+1);
132: mh_outstr((char *)cw->aux2,sfp->len+1,sfp);
133: }
134: /* deallocate the memory */
135: for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]);
136: mh_freeWorkArea();
137: return(cw);
138: }
139: /* Cumulative probability distribution function of the first eigenvalue of
140: Wishart matrix by Series */
141: struct cWishart *mh_cwishart_s(int m,int n,double beta[],double x0,
1.13 takayama 142: int approxDeg,double h, int dp, double x,
143: int automatic,double assigned_series_error,
144: int verbose)
145: {
1.6 takayama 146: int modep[]={1,0,0};
1.13 takayama 147: return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep,automatic,assigned_series_error,verbose));
1.1 takayama 148: }
149:
150: /* Cumulative probability distribution function of the first eigenvalue of
151: Wishart matrix by HGM */
152: struct cWishart *mh_cwishart_hgm(int m,int n,double beta[],double x0,
1.13 takayama 153: int approxDeg, double h, int dp , double x,
154: int automatic,double assigned_series_error,
155: int verbose)
1.1 takayama 156: {
1.6 takayama 157: int modep[]={1,1,0};
1.13 takayama 158: return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep,automatic,assigned_series_error,verbose));
1.1 takayama 159: }
160:
1.16 takayama 161: struct cWishart *mh_pFq_gen(int m,
162: int p, double a[],
163: int q, double b[],
164: int ef_type,
165: double beta[],double x0,
166: int approxDeg,double h, int dp, double x,int modep[],
167: int automatic,double assigned_series_error,
168: int verbose) {
169: /*
170: modep[0]. Do Koev-Edelman (ignored for now).
171: modep[1]. Do the HGM
172: modep[2]. Return modep[2] intermediate data.
173: */
174: struct SFILE *fp;
175: char swork[WSIZE];
176: char *argv[WSIZE];
177: int argc;
178: int i,rank;
179: char *comm;
180: struct MH_RESULT *rp;
181: struct SFILE *sfp;
182: struct cWishart *cw;
183: argv[0]="dummy";
184: argv[1] = "--bystring";
185: argv[2] = "--idata";
186: fp = mh_fopen("","w",0);
187: mh_fputs("%!version2.0\n",fp);
188: mh_fputs("%Mg=\n",fp);
189: sprintf(swork,"%d\n",m); mh_fputs(swork,fp);
190: rank = imypower(q+1,m);
191:
192: sprintf(swork,"%%p_pFq=%d ",p); mh_fputs(swork,fp);
193: for (i=0; i<p; i++) {
194: sprintf(swork,"%lg ",a[i]); mh_fputs(swork,fp);
195: }
196: mh_fputs("\n",fp);
197:
198: sprintf(swork,"%%q_pFq=%d ",q); mh_fputs(swork,fp);
199: for (i=0; i<q; i++) {
200: sprintf(swork,"%lg ",b[i]); mh_fputs(swork,fp);
201: }
202: mh_fputs("\n",fp);
203:
204: sprintf(swork,"%%ef_type=%d\n",ef_type); mh_fputs(swork,fp);
205:
206: mh_fputs("%Beta=\n",fp);
207: for (i=0; i<m; i++) {
208: sprintf(swork,"%lf\n",beta[i]); mh_fputs(swork,fp);
209: }
210: mh_fputs("%X0g=\n",fp); /* initial point */
211: sprintf(swork,"%lf\n",x0); mh_fputs(swork,fp);
212: if (h <= 0.0) {oxprintfe("h<=0.0, set to 0.1\n"); h=0.1;}
213: mh_fputs("%Hg=\n",fp);
214: sprintf(swork,"%lf\n",h); mh_fputs(swork,fp);
215: mh_fputs("%Dp=\n",fp);
216: sprintf(swork,"%d\n",dp); mh_fputs(swork,fp);
217: if (x <= x0) {oxprintfe("x <= x0, set to x=x0+10\n"); x=x0+10;}
218: mh_fputs("%Xng=\n",fp);
219: sprintf(swork,"%lf\n",x); mh_fputs(swork,fp);
220:
221: sprintf(swork,"%%automatic=%d\n",automatic); mh_fputs(swork,fp);
222: sprintf(swork,"%%assigned_series_error=%lg\n",assigned_series_error); mh_fputs(swork,fp);
223: sprintf(swork,"%%show_autosteps=%d\n",verbose); mh_fputs(swork,fp);
224:
225: comm = (char *)mh_malloc(fp->len +1);
226: mh_outstr(comm,fp->len+1,fp);
227: mh_fclose(fp);
228: argv[3] = comm;
1.17 takayama 229: if (MH_DEBUG2) oxprintf("comm=\n%s",comm);
1.16 takayama 230:
231: argv[4] = "--degree";
232: argv[5] = (char *)mh_malloc(128);
233: sprintf(argv[5],"%d",approxDeg);
234:
235: rp=jk_main(6,argv);
236: if (rp == NULL) {
237: oxprintfe("rp is NULL.\n"); return(NULL);
238: }
239: cw = new_cWishart(rank);
240: cw->x = rp->x;
241: cw->rank = rp->rank;
242: if (rank != cw->rank) {
243: oxprintfe("Rank error.\n"); return(NULL);
244: }
245: for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i];
246: sfp = (rp->sfpp)[0];
247: cw->aux = (char *) mh_malloc((sfp->len)+1);
248: mh_outstr((char *)(cw->aux),(sfp->len)+1,sfp);
249: /* todo, the following line seems to cause seg fault. */
250: /* deallocate the memory */
1.18 takayama 251: //for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]);
1.16 takayama 252: /* todo, mh_free_??(rp); free Iv's */
253: if (!modep[1]) return(cw);
254:
255: if (MH_DEBUG2) oxprintf("\n\n%s\n",(char *)cw->aux);
256: /* This output is strange. */
257: /* Starting HGM */
258: argv[3] = (char *)cw->aux;
259: argv[4] = "--dataf";
260: argv[5] = "dummy-dataf";
261: argc=6;
262: if (verbose) {
263: argv[6] = "--verbose"; ++argc;
264: }
265: rp = mh_main(argc,argv);
266: if (rp == NULL) {
267: oxprintfe("rp is NULL in the second step.\n"); return(NULL);
268: }
269: cw = new_cWishart(rank);
270: cw->x = rp->x;
271: cw->rank = rp->rank;
272: for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i];
273: sfp = (rp->sfpp)[0];
274: if (sfp) {
275: cw->aux = (char *) mh_malloc(sfp->len+1);
276: mh_outstr((char *)cw->aux,sfp->len+1,sfp);
277: }
278: sfp = (rp->sfpp)[1];
279: if (sfp) {
280: cw->aux2 = (char *) mh_malloc(sfp->len+1);
281: mh_outstr((char *)cw->aux2,sfp->len+1,sfp);
282: }
283: /* deallocate the memory */
284: for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]);
285: mh_freeWorkArea();
286: return(cw);
287: }
288:
289:
1.1 takayama 290: #ifdef STANDALONE
1.15 takayama 291: int main(int argc,char *argv[]) {
1.1 takayama 292: double beta[5]={1.0,2.0,3.0,4.0,5.0};
293: struct cWishart *cw;
1.7 takayama 294: struct SFILE *sfp;
295: char *s;
296: char str[1024];
297: double x;
1.12 takayama 298: int i,show;
1.13 takayama 299: int strategy=1;
1.12 takayama 300: double err[2]={-1.0,-1.0};
1.13 takayama 301: int verbose=0;
1.16 takayama 302: int testnew=0;
1.12 takayama 303: show=1;
304: for (i=1; i<argc; i++) {
305: if (strcmp(argv[i],"--strategy")==0) {
306: i++; sscanf(argv[i],"%d",&strategy);
307: }else if (strcmp(argv[i],"--abserr")==0) {
308: i++; sscanf(argv[i],"%lg",&(err[0]));
309: }else if (strcmp(argv[i],"--relerr")==0) {
310: i++; sscanf(argv[i],"%lg",&(err[1]));
311: }else if (strcmp(argv[i],"--quiet")==0) {
312: show=0;
1.13 takayama 313: }else if (strcmp(argv[i],"--verbose")==0) {
314: verbose=1;
1.16 takayama 315: }else if (strcmp(argv[i],"--testnew")==0) {
316: testnew=1;
1.12 takayama 317: }else{
1.16 takayama 318: oxprintfe("Unknown option.\n"); return(-1);
1.12 takayama 319: }
320: }
1.16 takayama 321: if (testnew) {
322: int tmodep[] = {1,1,0};
323: double a[] = {2.0,7.5};
324: double b[] = {4.5};
325: double beta[] = {1,2,4};
326: mh_set_strategy(strategy,err);
327: /* Testdata/new-2016-02-04-4-in.txt, Hashiguchi note 2016.02.03 */
328: cw=mh_pFq_gen(3, /* m */
329: 2,a,
330: 1,b,
331: 2, /* ef_type */
332: beta,0.1,
333: 8,0.001,1, 4,tmodep,
334: 1,1e-20,verbose);
335: if (cw != NULL) {
336: oxprintf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
337: /* oxprintf("%s",(char *)cw->aux); */
338: }
339: return(0);
340: }
1.12 takayama 341: mh_set_strategy(strategy,err);
1.13 takayama 342: cw=mh_cwishart_hgm(3,5,beta,0.3,7, 0.01,1,10, 1,1e-7,verbose);
1.1 takayama 343: if (cw != NULL) {
1.14 takayama 344: oxprintf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
345: /* oxprintf("%s",(char *)cw->aux); */
1.11 takayama 346: }
1.13 takayama 347: cw=mh_cwishart_hgm(4,5,beta,0.3,7, 0.01,1,10, 1,1e-7,verbose);
1.1 takayama 348: if (cw != NULL) {
1.14 takayama 349: oxprintf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
1.7 takayama 350: s = (char *)cw->aux;
1.14 takayama 351: /* oxprintf("%s",(char *)cw->aux); */
1.12 takayama 352: if (show) {
353: sfp = mh_fopen(s,"r",0);
354: while (mh_fgets(str,1024,sfp)) {
1.14 takayama 355: sscanf(str,"%lg",&x); oxprintf("%lg\n",x);
1.12 takayama 356: }
357: mh_fclose(sfp);
1.7 takayama 358: }
1.1 takayama 359: }
1.15 takayama 360: return(0);
1.1 takayama 361: }
1.15 takayama 362: int main1() {
1.1 takayama 363: double beta[5]={1.0,2.0,3.0,4.0,5.0};
364: struct cWishart *cw;
1.13 takayama 365: int verbose=1;
366: cw=mh_cwishart_s(3,5,beta,0.3,7, 0,0,0, 1,1e-7,verbose);
1.1 takayama 367: if (cw != NULL) {
1.14 takayama 368: oxprintf("%s",(char *)cw->aux);
1.1 takayama 369: }
1.13 takayama 370: cw=mh_cwishart_s(4,5,beta,0.3,7, 0,0,0, 1,1e-7,verbose);
1.1 takayama 371: if (cw != NULL) {
1.14 takayama 372: oxprintf("%s",(char *)cw->aux);
1.1 takayama 373: }
1.13 takayama 374: cw=mh_cwishart_s(5,5,beta,0.3,7, 0,0,0, 1,1e-7,verbose);
1.1 takayama 375: if (cw != NULL) {
1.14 takayama 376: oxprintf("%s",(char *)cw->aux);
1.1 takayama 377: }
1.15 takayama 378: return(0);
1.1 takayama 379: }
380: #endif
1.16 takayama 381:
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