version 1.8, 2013/03/05 07:03:37 |
version 1.16, 2016/02/13 05:55:09 |
|
|
/* $OpenXM: OpenXM/src/hgm/mh/src/mh.c,v 1.7 2013/03/05 06:35:54 takayama Exp $ */ |
/* $OpenXM: OpenXM/src/hgm/mh/src/mh.c,v 1.15 2015/04/02 01:11:13 takayama Exp $ */ |
#include <stdio.h> |
#include <stdio.h> |
|
#include <string.h> |
#include "sfile.h" |
#include "sfile.h" |
#include "mh.h" |
#include "mh.h" |
#define WSIZE 1024 |
#define WSIZE 1024 |
|
int MH_DEBUG2=0; |
extern int MH_DEBUG; |
extern int MH_DEBUG; |
static imypower(int x,int n) { |
extern int M_show_autosteps; |
|
static int imypower(int x,int n) { |
int a,i; |
int a,i; |
a = 1; |
a = 1; |
for (i=0; i<n; i++) a = a*x; |
for (i=0; i<n; i++) a = a*x; |
Line 22 struct cWishart *new_cWishart(int rank) { |
|
Line 25 struct cWishart *new_cWishart(int rank) { |
|
} |
} |
|
|
struct cWishart *mh_cwishart_gen(int m,int n,double beta[],double x0, |
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[0]. Do Koev-Edelman (ignored for now). |
modep[1]. Do the HGM |
modep[1]. Do the HGM |
Line 31 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
Line 37 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
struct SFILE *fp; |
struct SFILE *fp; |
char swork[WSIZE]; |
char swork[WSIZE]; |
char *argv[WSIZE]; |
char *argv[WSIZE]; |
|
int argc; |
int i,rank; |
int i,rank; |
char *comm; |
char *comm; |
struct MH_RESULT *rp; |
struct MH_RESULT *rp; |
Line 54 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
Line 61 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
mh_fputs("%%Iv\n",fp); /* initial values, dummy */ |
mh_fputs("%%Iv\n",fp); /* initial values, dummy */ |
for (i=0; i<rank; i++) mh_fputs("0.0\n",fp); |
for (i=0; i<rank; i++) mh_fputs("0.0\n",fp); |
mh_fputs("%%Ef=\n1.0\n",fp); /* Below are dummy values */ |
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); |
mh_fputs("%%Hg=\n",fp); |
sprintf(swork,"%lf\n",h); mh_fputs(swork,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); |
mh_fputs("%%Dp=\n",fp); |
sprintf(swork,"%d\n",dp); mh_fputs(swork,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); |
mh_fputs("%%Xng=\n",fp); |
sprintf(swork,"%lf\n",x); mh_fputs(swork,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); |
comm = (char *)mh_malloc(fp->len +1); |
mh_outstr(comm,fp->len+1,fp); |
mh_outstr(comm,fp->len+1,fp); |
mh_fclose(fp); |
mh_fclose(fp); |
Line 75 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
Line 86 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
|
|
rp=jk_main(6,argv); |
rp=jk_main(6,argv); |
if (rp == NULL) { |
if (rp == NULL) { |
fprintf(stderr,"rp is NULL.\n"); return(NULL); |
oxprintfe("rp is NULL.\n"); return(NULL); |
} |
} |
cw = new_cWishart(rank); |
cw = new_cWishart(rank); |
cw->x = rp->x; |
cw->x = rp->x; |
cw->rank = rp->rank; |
cw->rank = rp->rank; |
if (rank != cw->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]; |
for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i]; |
sfp = (rp->sfpp)[0]; |
sfp = (rp->sfpp)[0]; |
Line 93 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
Line 104 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
/* todo, mh_free_??(rp); free Iv's */ |
/* todo, mh_free_??(rp); free Iv's */ |
if (!modep[1]) return(cw); |
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. */ |
/* This output is strange. */ |
/* Starting HGM */ |
/* Starting HGM */ |
argv[3] = (char *)cw->aux; |
argv[3] = (char *)cw->aux; |
argv[4] = "--dataf"; |
argv[4] = "--dataf"; |
argv[5] = "dummy-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) { |
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 = new_cWishart(rank); |
cw->x = rp->x; |
cw->x = rp->x; |
Line 125 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
Line 140 struct cWishart *mh_cwishart_gen(int m,int n,double be |
|
/* Cumulative probability distribution function of the first eigenvalue of |
/* Cumulative probability distribution function of the first eigenvalue of |
Wishart matrix by Series */ |
Wishart matrix by Series */ |
struct cWishart *mh_cwishart_s(int m,int n,double beta[],double x0, |
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}; |
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 |
/* Cumulative probability distribution function of the first eigenvalue of |
Wishart matrix by HGM */ |
Wishart matrix by HGM */ |
struct cWishart *mh_cwishart_hgm(int m,int n,double beta[],double x0, |
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}; |
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)); |
} |
} |
|
|
|
struct cWishart *mh_pFq_gen(int m, |
|
int p, double a[], |
|
int q, double b[], |
|
int ef_type, |
|
double beta[],double x0, |
|
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 |
|
modep[2]. Return modep[2] intermediate data. |
|
*/ |
|
struct SFILE *fp; |
|
char swork[WSIZE]; |
|
char *argv[WSIZE]; |
|
int argc; |
|
int i,rank; |
|
char *comm; |
|
struct MH_RESULT *rp; |
|
struct SFILE *sfp; |
|
struct cWishart *cw; |
|
argv[0]="dummy"; |
|
argv[1] = "--bystring"; |
|
argv[2] = "--idata"; |
|
fp = mh_fopen("","w",0); |
|
mh_fputs("%!version2.0\n",fp); |
|
mh_fputs("%Mg=\n",fp); |
|
sprintf(swork,"%d\n",m); mh_fputs(swork,fp); |
|
rank = imypower(q+1,m); |
|
|
|
sprintf(swork,"%%p_pFq=%d ",p); mh_fputs(swork,fp); |
|
for (i=0; i<p; i++) { |
|
sprintf(swork,"%lg ",a[i]); mh_fputs(swork,fp); |
|
} |
|
mh_fputs("\n",fp); |
|
|
|
sprintf(swork,"%%q_pFq=%d ",q); mh_fputs(swork,fp); |
|
for (i=0; i<q; i++) { |
|
sprintf(swork,"%lg ",b[i]); mh_fputs(swork,fp); |
|
} |
|
mh_fputs("\n",fp); |
|
|
|
sprintf(swork,"%%ef_type=%d\n",ef_type); mh_fputs(swork,fp); |
|
|
|
mh_fputs("%Beta=\n",fp); |
|
for (i=0; i<m; i++) { |
|
sprintf(swork,"%lf\n",beta[i]); mh_fputs(swork,fp); |
|
} |
|
mh_fputs("%X0g=\n",fp); /* initial point */ |
|
sprintf(swork,"%lf\n",x0); mh_fputs(swork,fp); |
|
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) {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) {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); |
|
argv[3] = comm; |
|
if (MH_DEBUG2) printf("comm=\n%s",comm); |
|
|
|
argv[4] = "--degree"; |
|
argv[5] = (char *)mh_malloc(128); |
|
sprintf(argv[5],"%d",approxDeg); |
|
|
|
rp=jk_main(6,argv); |
|
if (rp == NULL) { |
|
oxprintfe("rp is NULL.\n"); return(NULL); |
|
} |
|
cw = new_cWishart(rank); |
|
cw->x = rp->x; |
|
cw->rank = rp->rank; |
|
if (rank != cw->rank) { |
|
oxprintfe("Rank error.\n"); return(NULL); |
|
} |
|
for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i]; |
|
sfp = (rp->sfpp)[0]; |
|
cw->aux = (char *) mh_malloc((sfp->len)+1); |
|
mh_outstr((char *)(cw->aux),(sfp->len)+1,sfp); |
|
/* todo, the following line seems to cause seg fault. */ |
|
/* deallocate the memory */ |
|
for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]); |
|
/* todo, mh_free_??(rp); free Iv's */ |
|
if (!modep[1]) return(cw); |
|
|
|
if (MH_DEBUG2) 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"; |
|
argc=6; |
|
if (verbose) { |
|
argv[6] = "--verbose"; ++argc; |
|
} |
|
rp = mh_main(argc,argv); |
|
if (rp == NULL) { |
|
oxprintfe("rp is NULL in the second step.\n"); return(NULL); |
|
} |
|
cw = new_cWishart(rank); |
|
cw->x = rp->x; |
|
cw->rank = rp->rank; |
|
for (i=0; i<cw->rank; i++) (cw->f)[i] = (rp->y)[i]; |
|
sfp = (rp->sfpp)[0]; |
|
if (sfp) { |
|
cw->aux = (char *) mh_malloc(sfp->len+1); |
|
mh_outstr((char *)cw->aux,sfp->len+1,sfp); |
|
} |
|
sfp = (rp->sfpp)[1]; |
|
if (sfp) { |
|
cw->aux2 = (char *) mh_malloc(sfp->len+1); |
|
mh_outstr((char *)cw->aux2,sfp->len+1,sfp); |
|
} |
|
/* deallocate the memory */ |
|
for (i=0; i<rp->size; i++) mh_fclose((rp->sfpp)[i]); |
|
mh_freeWorkArea(); |
|
return(cw); |
|
} |
|
|
|
|
#ifdef STANDALONE |
#ifdef STANDALONE |
main() { |
int main(int argc,char *argv[]) { |
double beta[5]={1.0,2.0,3.0,4.0,5.0}; |
double beta[5]={1.0,2.0,3.0,4.0,5.0}; |
struct cWishart *cw; |
struct cWishart *cw; |
struct SFILE *sfp; |
struct SFILE *sfp; |
char *s; |
char *s; |
char str[1024]; |
char str[1024]; |
double x; |
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; |
|
int testnew=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 if (strcmp(argv[i],"--testnew")==0) { |
|
testnew=1; |
|
}else{ |
|
oxprintfe("Unknown option.\n"); return(-1); |
|
} |
|
} |
|
if (testnew) { |
|
int tmodep[] = {1,1,0}; |
|
double a[] = {2.0,7.5}; |
|
double b[] = {4.5}; |
|
double beta[] = {1,2,4}; |
|
mh_set_strategy(strategy,err); |
|
/* Testdata/new-2016-02-04-4-in.txt, Hashiguchi note 2016.02.03 */ |
|
cw=mh_pFq_gen(3, /* m */ |
|
2,a, |
|
1,b, |
|
2, /* ef_type */ |
|
beta,0.1, |
|
8,0.001,1, 4,tmodep, |
|
1,1e-20,verbose); |
|
if (cw != NULL) { |
|
oxprintf("x=%lf, y=%lf\n",cw->x,(cw->f)[0]); |
|
/* oxprintf("%s",(char *)cw->aux); */ |
|
} |
|
return(0); |
|
} |
|
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) { |
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]); |
/* printf("%s",(char *)cw->aux); */ |
/* 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) { |
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; |
s = (char *)cw->aux; |
/* printf("%s",(char *)cw->aux); */ |
/* oxprintf("%s",(char *)cw->aux); */ |
sfp = mh_fopen(s,"r",0); |
if (show) { |
while (mh_fgets(str,1024,sfp)) { |
sfp = mh_fopen(s,"r",0); |
sscanf(str,"%lg",&x); printf("%lg\n",x); |
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}; |
double beta[5]={1.0,2.0,3.0,4.0,5.0}; |
struct cWishart *cw; |
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) { |
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) { |
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) { |
if (cw != NULL) { |
printf("%s",(char *)cw->aux); |
oxprintf("%s",(char *)cw->aux); |
} |
} |
|
return(0); |
} |
} |
#endif |
#endif |
|
|