OpenXM/src/hgm/mh/src/wmain.c - diff

Return to wmain.c CVS log

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

Diff for /OpenXM/src/hgm/mh/src/wmain.c between version 1.21 and 1.22

version 1.21, 2014/03/20 10:58:37

version 1.22, 2015/03/24 05:59:43

Line 1

$OpenXM: OpenXM/src/hgm/mh/src/wmain.c,v 1.20 2014/03/15 00:43:47 takayama Exp $

$OpenXM: OpenXM/src/hgm/mh/src/wmain.c,v 1.21 2014/03/20 10:58:37 takayama Exp $

License: LGPL

#include <stdio.h>

Line 9

#include "sfile.h"

#include "mh.h"

#define SMAX 4096

#define inci(i) { i++; if (i >= argc) { fprintf(stderr,"Option argument is not given.\n"); return(NULL); }}

#define inci(i) { i++; if (i >= argc) { oxprintfe("Option argument is not given.\n"); return(NULL); }}

int MH_deallocate=0;

Line 147 struct MH_RESULT *mh_main(int argc,char *argv[]) {

}else if (strcmp(argv[i],"--relerr")==0) {

i++; sscanf(argv[i],"%lg",&MH_relerr);

}else {

fprintf(stderr,"Unknown option %s\n",argv[i]);

oxprintfe("Unknown option %s\n",argv[i]);

mh_usage();

return(rp);

}

Line 164 struct MH_RESULT *mh_main(int argc,char *argv[]) {

if (MH_strategy) {

if (MH_abserr > SIGDIGIT_DEFAULT*myabs(y0[0])) {

MH_success = 0;

fprintf(stderr,"%%%%Warning, abserr seems not to be small enough, abserr=%lg, y[0]=%lg\n",MH_abserr,y0[0]);

oxprintfe("%%%%Warning, abserr seems not to be small enough, abserr=%lg, y[0]=%lg\n",MH_abserr,y0[0]);

}else{

MH_success = 1;

}

Line 173 struct MH_RESULT *mh_main(int argc,char *argv[]) {

}

MH_estimated_start_step = estimateHg(MH_Mg,MH_Beta,MH_X0g);

if (MH_Verbose) showParam();

if (MH_Verbose) {for (i=0; i<rank; i++) printf("%lf\n",y0[i]); }

if (MH_Verbose) {for (i=0; i<rank; i++) oxprintf("%lf\n",y0[i]); }

*rp=mh_rkmain(x0,y0,xn);

return(rp);

}

int mh_usage() {

fprintf(stderr,"Usages:\n");

oxprintfe("Usages:\n");

fprintf(stderr,"hgm_w-n [--idata input_data_file --gnuplotf gnuplot_file_name\n");

oxprintfe("hgm_w-n [--idata input_data_file --gnuplotf gnuplot_file_name\n");

fprintf(stderr," --dataf output_data_file --raw --xmax xmax --test m --step h]\n");

oxprintfe(" --dataf output_data_file --raw --xmax xmax --test m --step h]\n");

fprintf(stderr,"[ --95 --verbose] \n");

oxprintfe("[ --95 --verbose] \n");

fprintf(stderr,"[ --strategy s --abserr ae --relerr re] \n");

oxprintfe("[ --strategy s --abserr ae --relerr re] \n");

fprintf(stderr,"s:0 rk, s:1 adaptive, s:2 adaptive&multiply, see rk.c for the default value of ae and re.\n");

oxprintfe("s:0 rk, s:1 adaptive, s:2 adaptive&multiply, see rk.c for the default value of ae and re.\n");

fprintf(stderr,"strategy default = %d\n",MH_strategy);

oxprintfe("strategy default = %d\n",MH_strategy);

fprintf(stderr,"--raw does not add data parameters to the output_data_file.\n");

oxprintfe("--raw does not add data parameters to the output_data_file.\n");

fprintf(stderr,"\nThe command hgm_w-n [options] evaluates Pr({y | y<xmax}), which is the cumulative distribution function of the largest root of the m by m Wishart matrix with n degrees of freedom and the covariantce matrix sigma.\n");

oxprintfe("\nThe command hgm_w-n [options] evaluates Pr({y | y<xmax}), which is the cumulative distribution function of the largest root of the m by m Wishart matrix with n degrees of freedom and the covariantce matrix sigma.\n");

fprintf(stderr,"All the eigenvalues of sigma must be simple.\n");

oxprintfe("All the eigenvalues of sigma must be simple.\n");

fprintf(stderr,"Parameters are specified by the input_data_file.\n");

oxprintfe("Parameters are specified by the input_data_file.\n");

fprintf(stderr,"Parameters are redefined when they appear more than once in the idata file and the command line options.\n");

oxprintfe("Parameters are redefined when they appear more than once in the idata file and the command line options.\n");

fprintf(stderr,"The format of the input_data_file, which should be generated by the command hgm_jack-n.\n");

oxprintfe("The format of the input_data_file, which should be generated by the command hgm_jack-n.\n");

fprintf(stderr," MH_Mg: m, MH_Beta: beta=sigma^(-1)/2 (diagonized), MH_Ng: n, MH_X0g: starting value of x,\n");

oxprintfe(" MH_Mg: m, MH_Beta: beta=sigma^(-1)/2 (diagonized), MH_Ng: n, MH_X0g: starting value of x,\n");

fprintf(stderr," Iv: initial values at MH_X0g*MH_Beta (see our paper how to order them), \n");

oxprintfe(" Iv: initial values at MH_X0g*MH_Beta (see our paper how to order them), \n");

fprintf(stderr," Ef: a scalar factor to the initial value. It may set to 1.\n");

oxprintfe(" Ef: a scalar factor to the initial value. It may set to 1.\n");

fprintf(stderr," MH_Hg: h (step size),\n");

oxprintfe(" MH_Hg: h (step size),\n");

fprintf(stderr," MH_Dp: output data is stored in every MH_Dp steps when output_data_file is specified.\n");

oxprintfe(" MH_Dp: output data is stored in every MH_Dp steps when output_data_file is specified.\n");

fprintf(stderr," Xng: terminating value of x.\n\n");

oxprintfe(" Xng: terminating value of x.\n\n");

fprintf(stderr,"--95: output the 95%% point. --verbose: verbose mode.\n");

oxprintfe("--95: output the 95%% point. --verbose: verbose mode.\n");

fprintf(stderr,"The line started with %%%% or # is a comment line.\n");

oxprintfe("The line started with %%%% or # is a comment line.\n");

fprintf(stderr,"An example format of the input_data_file can be obtained by executing hgm_jack-n with no option.\n");

oxprintfe("An example format of the input_data_file can be obtained by executing hgm_jack-n with no option.\n");

fprintf(stderr,"When --idata option is used, this command is quiet. Use --verbose option if you want to see some messages.\n");

oxprintfe("When --idata option is used, this command is quiet. Use --verbose option if you want to see some messages.\n");

fprintf(stderr,"\nExamples:\n");

oxprintfe("\nExamples:\n");

fprintf(stderr,"[1] ./hgm_w-n \n");

oxprintfe("[1] ./hgm_w-n \n");

fprintf(stderr,"[2] ./hgm_w-n --xmax 20\n");

oxprintfe("[2] ./hgm_w-n --xmax 20\n");

fprintf(stderr,"[3] ./hgm_w-n --test 6\n");

oxprintfe("[3] ./hgm_w-n --test 6\n");

fprintf(stderr," A test run in Mg=6.\n");

oxprintfe(" A test run in Mg=6.\n");

fprintf(stderr,"[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt >t.txt\n");

oxprintfe("[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt >t.txt\n");

fprintf(stderr," ./hgm_w-n --idata t.txt --gnuplotf test-g --verbose\n");

oxprintfe(" ./hgm_w-n --idata t.txt --gnuplotf test-g --verbose\n");

fprintf(stderr," gnuplot -persist <test-g-gp.txt\n");

oxprintfe(" gnuplot -persist <test-g-gp.txt\n");

fprintf(stderr," tmp-idata3.txt is a sample input data distributed with this file.\n");

oxprintfe(" tmp-idata3.txt is a sample input data distributed with this file.\n");

fprintf(stderr," test-g-gp.txt is an input file of the gnuplot\n");

oxprintfe(" test-g-gp.txt is an input file of the gnuplot\n");

fprintf(stderr," test-g is the table of x and the values of Pr({y | y<x}).\n");

oxprintfe(" test-g is the table of x and the values of Pr({y | y<x}).\n");

}

static int setParamTest() {

Line 259 static int next(struct SFILE *sfp,char *s,char *msg) {

s[0] = '%';

while (s[0] == '%') {

if (!mh_fgets(s,SMAX,sfp)) {

fprintf(stderr,"Data format error at %s\n",msg);

oxprintfe("Data format error at %s\n",msg);

mh_exit(-1);

}

if (s[0] != '%') return(0);

Line 282 static int setParam(char *fname) {

if (fname == NULL) return(setParamDefault());

if ((fp=mh_fopen(fname,"r",MH_byFile)) == NULL) {

fprintf(stderr,"File %s is not found.\n",fname);

oxprintfe("File %s is not found.\n",fname);

mh_exit(-1);

}

next(fp,s,"MH_Mg(m)");

Line 324 static int setParam(char *fname) {

while ((tk = mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_EOF) {

/* expect ID */

if (tk.type != MH_TOKEN_ID) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

if ((strcmp(s,"abserr")==0) || (strcmp(s,"abserror")==0)) {

if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

MH_abserr = tk.dval;

continue;

}

if ((strcmp(s,"relerr")==0) || (strcmp(s,"relerror")==0)) {

if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_DOUBLE) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

MH_relerr = tk.dval;

continue;

}

if (strcmp(s,"strategy")==0) {

if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {

fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);

oxprintfe("Syntax error at %s\n",s); mh_exit(-1);

}

MH_strategy = tk.ival;

continue;

}

fprintf(stderr,"Unknown ID at %s\n",s); mh_exit(-1);

oxprintfe("Unknown ID at %s\n",s); mh_exit(-1);

}

mh_fclose(fp);

Line 368 static int showParam() {

extern double MH_abserr;

extern double MH_relerr;

rank = mypower(2,MH_Mg);

printf("%%MH_Mg=\n%d\n",MH_Mg);

oxprintf("%%MH_Mg=\n%d\n",MH_Mg);

for (i=0; i<MH_Mg; i++) {

printf("%%MH_Beta[%d]=\n%lf\n",i,MH_Beta[i]);

oxprintf("%%MH_Beta[%d]=\n%lf\n",i,MH_Beta[i]);

}

printf("%%MH_Ng=\n%lf\n",*MH_Ng);

oxprintf("%%MH_Ng=\n%lf\n",*MH_Ng);

printf("%%MH_X0g=\n%lf\n",MH_X0g);

oxprintf("%%MH_X0g=\n%lf\n",MH_X0g);

for (i=0; i<rank; i++) {

printf("%%Iv[%d]=\n%lg\n",i,Iv[i]);

oxprintf("%%Iv[%d]=\n%lg\n",i,Iv[i]);

}

printf("%%Ef=\n%lf\n",Ef);

oxprintf("%%Ef=\n%lf\n",Ef);

printf("%%MH_Hg=\n%lf\n",MH_Hg);

oxprintf("%%MH_Hg=\n%lf\n",MH_Hg);

printf("%%MH_Dp=\n%d\n",MH_Dp);

oxprintf("%%MH_Dp=\n%d\n",MH_Dp);

printf("%%Xng=\n%lf\n",Xng);

oxprintf("%%Xng=\n%lf\n",Xng);

printf("%%strategy=%d\n",MH_strategy);

oxprintf("%%strategy=%d\n",MH_strategy);

printf("%%abserr=%lg, %%relerr=%lg\n",MH_abserr,MH_relerr);

oxprintf("%%abserr=%lg, %%relerr=%lg\n",MH_abserr,MH_relerr);

printf("#MH_success=%d\n",MH_success);

oxprintf("#MH_success=%d\n",MH_success);

printf("#MH_coeff_max=%lg\n",MH_coeff_max);

oxprintf("#MH_coeff_max=%lg\n",MH_coeff_max);

printf("#MH_estimated_start_step=%lg\n",MH_estimated_start_step);

oxprintf("#MH_estimated_start_step=%lg\n",MH_estimated_start_step);

}

static double estimateHg(int m, double beta[],double x0) {

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>