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