Return to mh.c CVS log

Up to [local] / OpenXM / src / hgm / mh / src

Annotation of OpenXM/src/hgm/mh/src/mh.c, Revision 1.3

1.3     ! takayama    1: /* $OpenXM: OpenXM/src/hgm/mh/src/mh.c,v 1.2 2013/02/23 06:01:53 takayama Exp $ */
1.1       takayama    2: #include <stdio.h>
                      3: #include "sfile.h"
                      4: #include "mh.h"
                      5: #define WSIZE 1024
                      6:
                      7: static imypower(int x,int n) {
                      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.1       takayama   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) {
1.3     ! takayama  119:   return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,0));
1.1       takayama  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: {
1.3     ! takayama  127:   return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,1));
1.1       takayama  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