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.2

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