===================================================================
RCS file: /home/cvs/OpenXM/src/hgm/mh/src/mh.c,v
retrieving revision 1.11
retrieving revision 1.15
diff -u -p -r1.11 -r1.15
--- OpenXM/src/hgm/mh/src/mh.c	2013/03/07 07:02:18	1.11
+++ OpenXM/src/hgm/mh/src/mh.c	2015/04/02 01:11:13	1.15
@@ -1,9 +1,11 @@
-/* $OpenXM: OpenXM/src/hgm/mh/src/mh.c,v 1.10 2013/03/07 05:23:31 takayama Exp $ */
+/* $OpenXM: OpenXM/src/hgm/mh/src/mh.c,v 1.14 2015/03/24 05:59:43 takayama Exp $ */
 #include <stdio.h>
+#include <string.h>
 #include "sfile.h"
 #include "mh.h"
 #define WSIZE 1024
 extern int MH_DEBUG;
+extern int M_show_autosteps;
 static int imypower(int x,int n) {
   int a,i;
   a = 1;
@@ -22,7 +24,10 @@ struct cWishart *new_cWishart(int rank) {
 }
 
 struct cWishart *mh_cwishart_gen(int m,int n,double beta[],double x0,
-                              int approxDeg,double h, int dp, double x,int modep[]) {
+				 int approxDeg,double h, int dp, double x,int modep[],
+                               int automatic,double assigned_series_error,
+                               int verbose)
+{
   /*
      modep[0]. Do Koev-Edelman (ignored for now).
      modep[1]. Do the HGM
@@ -31,6 +36,7 @@ struct cWishart *mh_cwishart_gen(int m,int n,double be
   struct SFILE *fp;
   char swork[WSIZE];
   char *argv[WSIZE];
+  int argc;
   int i,rank;
   char *comm;
   struct MH_RESULT *rp;
@@ -54,16 +60,20 @@ struct cWishart *mh_cwishart_gen(int m,int n,double be
   mh_fputs("%%Iv\n",fp); /* initial values, dummy */
   for (i=0; i<rank; i++) mh_fputs("0.0\n",fp);
   mh_fputs("%%Ef=\n1.0\n",fp); /* Below are dummy values */
-  if (h <= 0.0) {fprintf(stderr,"h<=0.0, set to 0.1\n"); h=0.1;}
+  if (h <= 0.0) {oxprintfe("h<=0.0, set to 0.1\n"); h=0.1;}
   mh_fputs("%%Hg=\n",fp); 
   sprintf(swork,"%lf\n",h); mh_fputs(swork,fp);
-  if (dp < 1) {fprintf(stderr,"dp<1, set to 1\n"); dp=1;}
+  if (dp < 1) {oxprintfe("dp<1, set to 1\n"); dp=1;}
   mh_fputs("%%Dp=\n",fp); 
   sprintf(swork,"%d\n",dp); mh_fputs(swork,fp);
-  if (x <= x0) {fprintf(stderr,"x <= x0, set to x=x0+10\n"); x=x0+10;}
+  if (x <= x0) {oxprintfe("x <= x0, set to x=x0+10\n"); x=x0+10;}
   mh_fputs("%%Xng=\n",fp);
   sprintf(swork,"%lf\n",x); mh_fputs(swork,fp);
 
+  sprintf(swork,"%%automatic=%d\n",automatic); mh_fputs(swork,fp);
+  sprintf(swork,"%%assigned_series_error=%lg\n",assigned_series_error); mh_fputs(swork,fp);
+  sprintf(swork,"%%show_autosteps=%d\n",verbose); mh_fputs(swork,fp);
+
   comm = (char *)mh_malloc(fp->len +1);
   mh_outstr(comm,fp->len+1,fp);
   mh_fclose(fp);
@@ -75,13 +85,13 @@ struct cWishart *mh_cwishart_gen(int m,int n,double be
 
   rp=jk_main(6,argv); 
   if (rp == NULL) {
-    fprintf(stderr,"rp is NULL.\n"); return(NULL);
+    oxprintfe("rp is NULL.\n"); return(NULL);
   }
   cw = new_cWishart(rank);
   cw->x = rp->x;
   cw->rank = rp->rank;
   if (rank !=  cw->rank) {
-    fprintf(stderr,"Rank error.\n"); return(NULL);
+    oxprintfe("Rank error.\n"); return(NULL);
   }
   for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i];  
   sfp = (rp->sfpp)[0]; 
@@ -93,15 +103,19 @@ struct cWishart *mh_cwishart_gen(int m,int n,double be
   /* todo, mh_free_??(rp);  free Iv's */
   if (!modep[1]) return(cw);
 
-  if (MH_DEBUG) printf("\n\n%s\n",(char *)cw->aux);
+  if (MH_DEBUG) oxprintf("\n\n%s\n",(char *)cw->aux);
   /* This output is strange. */
   /* Starting HGM */
   argv[3] = (char *)cw->aux;
   argv[4] = "--dataf";
   argv[5] = "dummy-dataf";
-  rp = mh_main(6,argv);
+  argc=6;
+  if (verbose) {
+    argv[6] = "--verbose"; ++argc;
+  }
+  rp = mh_main(argc,argv);
   if (rp == NULL) {
-    fprintf(stderr,"rp is NULL in the second step.\n"); return(NULL);
+    oxprintfe("rp is NULL in the second step.\n"); return(NULL);
   }
   cw = new_cWishart(rank);
   cw->x = rp->x;
@@ -125,59 +139,90 @@ struct cWishart *mh_cwishart_gen(int m,int n,double be
 /* Cumulative probability distribution function of the first eigenvalue of
    Wishart matrix by Series */
 struct cWishart *mh_cwishart_s(int m,int n,double beta[],double x0,
-                               int approxDeg,double h, int dp, double x) {
+                               int approxDeg,double h, int dp, double x,
+                               int automatic,double assigned_series_error,
+                               int verbose)
+{
   int modep[]={1,0,0};
-  return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep));
+  return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep,automatic,assigned_series_error,verbose));
 }
 
 /* Cumulative probability distribution function of the first eigenvalue of
    Wishart matrix by HGM */
 struct cWishart *mh_cwishart_hgm(int m,int n,double beta[],double x0,
-                                 int approxDeg, double h, int dp , double x)
+                                 int approxDeg, double h, int dp , double x,
+                               int automatic,double assigned_series_error,
+                               int verbose)
 {
   int modep[]={1,1,0};
-  return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep));
+  return(mh_cwishart_gen(m,n,beta,x0,approxDeg,h,dp,x,modep,automatic,assigned_series_error,verbose));
 }
 
 #ifdef STANDALONE
-main() {
+int main(int argc,char *argv[]) {
   double beta[5]={1.0,2.0,3.0,4.0,5.0};
   struct cWishart *cw;
   struct SFILE *sfp;
   char *s;
   char str[1024];
   double x;
-  cw=mh_cwishart_hgm(3,5,beta,0.3,7,  0.01,1,10);
+  int i,show;
+  int strategy=1;
+  double err[2]={-1.0,-1.0};
+  int verbose=0;
+  show=1;
+  for (i=1; i<argc; i++) {
+	if (strcmp(argv[i],"--strategy")==0) {
+	  i++; sscanf(argv[i],"%d",&strategy);
+	}else if (strcmp(argv[i],"--abserr")==0) {
+	  i++; sscanf(argv[i],"%lg",&(err[0]));
+	}else if (strcmp(argv[i],"--relerr")==0) {
+	  i++; sscanf(argv[i],"%lg",&(err[1]));
+	}else if (strcmp(argv[i],"--quiet")==0) {
+	  show=0;
+	}else if (strcmp(argv[i],"--verbose")==0) {
+	  verbose=1;
+	}else{
+	  oxprintfe("Unknown option.\n");
+	}
+  }
+  mh_set_strategy(strategy,err);
+  cw=mh_cwishart_hgm(3,5,beta,0.3,7,  0.01,1,10, 1,1e-7,verbose);
   if (cw != NULL) {
-    printf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
-    /* printf("%s",(char *)cw->aux); */
+    oxprintf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
+    /* oxprintf("%s",(char *)cw->aux); */
   }
-  cw=mh_cwishart_hgm(4,5,beta,0.3,7,  0.01,1,10);
+  cw=mh_cwishart_hgm(4,5,beta,0.3,7,  0.01,1,10, 1,1e-7,verbose);
   if (cw != NULL) {
-    printf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
+    oxprintf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]);
     s = (char *)cw->aux;
-    /* printf("%s",(char *)cw->aux); */
-    sfp = mh_fopen(s,"r",0);
-    while (mh_fgets(str,1024,sfp)) {
-      sscanf(str,"%lg",&x); printf("%lg\n",x);
+    /* oxprintf("%s",(char *)cw->aux); */
+    if (show) {
+      sfp = mh_fopen(s,"r",0);
+      while (mh_fgets(str,1024,sfp)) {
+        sscanf(str,"%lg",&x); oxprintf("%lg\n",x);
+      }
+      mh_fclose(sfp);
     }
-    mh_fclose(sfp);
   }
+  return(0);
 }
-main1() {
+int main1() {
   double beta[5]={1.0,2.0,3.0,4.0,5.0};
   struct cWishart *cw;
-  cw=mh_cwishart_s(3,5,beta,0.3,7,  0,0,0);
+  int verbose=1;
+  cw=mh_cwishart_s(3,5,beta,0.3,7,  0,0,0, 1,1e-7,verbose);
   if (cw != NULL) {
-    printf("%s",(char *)cw->aux);
+    oxprintf("%s",(char *)cw->aux);
   }
-  cw=mh_cwishart_s(4,5,beta,0.3,7,  0,0,0);
+  cw=mh_cwishart_s(4,5,beta,0.3,7,  0,0,0, 1,1e-7,verbose);
   if (cw != NULL) {
-    printf("%s",(char *)cw->aux);
+    oxprintf("%s",(char *)cw->aux);
   }
-  cw=mh_cwishart_s(5,5,beta,0.3,7,  0,0,0);
+  cw=mh_cwishart_s(5,5,beta,0.3,7,  0,0,0, 1,1e-7,verbose);
   if (cw != NULL) {
-    printf("%s",(char *)cw->aux);
+    oxprintf("%s",(char *)cw->aux);
   }
+  return(0);
 }
 #endif