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