Return to rk.c CVS log

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

Annotation of OpenXM/src/hgm/mh/src/rk.c, Revision 1.16

1.1       takayama    1: /*
1.7       takayama    2:   License: LGPL
                      3:   Ref: Copied from this11/misc-2011/A1/wishart/Prog
1.16    ! takayama    4:   $OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.15 2016/02/13 02:18:59 takayama Exp $
1.7       takayama    5: */
1.1       takayama    6: #include <stdio.h>
                      7: #include <stdlib.h>
                      8: #include <math.h>
                      9: #include <string.h>
1.4       takayama   10: #include "sfile.h"
1.8       takayama   11: #include "mh.h"
1.1       takayama   12:
1.9       takayama   13: #include "t-gsl_errno.h"
                     14: #include "t-gsl_odeiv.h"
1.1       takayama   15: extern int MH_RANK;
                     16: extern int MH_M;
1.15      takayama   17: extern MH_RF mh_rf;
1.1       takayama   18:
                     19: char *MH_Gfname;
                     20: char *MH_Dfname;
                     21: extern int MH_Mg;
1.2       takayama   22: extern double *MH_Beta; extern double *MH_Ng; extern double MH_X0g;
                     23: extern int MH_RawName;
                     24: double MH_Hg = 0.001;
1.4       takayama   25: static struct SFILE *Gf =NULL;
                     26: static struct SFILE *Df =NULL;
1.2       takayama   27: int MH_P95=0;
                     28: int MH_Verbose=0;
1.11      takayama   29: int MH_strategy=STRATEGY_DEFAULT;  /* 0: rk4, 1:a-rk4 (adaptive) */
                     30: double MH_abserr = 1e-10;
1.12      takayama   31: double MH_relerr = MH_RELERR_DEFAULT ;
1.11      takayama   32: /* D_i = MH_abserr + MH_relerr*|y_i|
                     33:    If observed error E > D*110/100, the stepsize will be decreased.
                     34:    cf. GSL, 26.3 Adaptive Step-size control.
                     35: */
1.9       takayama   36:
                     37: int mh_rf_for_gsl(double t,const double y[],double f[],void *params) {
1.15      takayama   38:   (*mh_rf)(t,(double *)y,MH_RANK,f,MH_RANK);
1.9       takayama   39:   return(GSL_SUCCESS);
                     40: }
                     41:
1.6       takayama   42: static int mypower(int x,int n) {
1.1       takayama   43:   int a,i;
                     44:   a = 1;
                     45:   for (i=0; i<n; i++) a = a*x;
                     46:   return(a);
                     47: }
1.6       takayama   48: int mh_gopen_file() {
1.1       takayama   49:   FILE *fp;
                     50:   char fname[1024];
                     51:   int i;
1.4       takayama   52:   extern int MH_byFile;
1.1       takayama   53:   Gf=NULL;
                     54:   Df=NULL;
1.2       takayama   55:   if (MH_Verbose) {
1.13      takayama   56:     oxprintfe("MH_Gfname=%s\n",MH_Gfname);
                     57:     oxprintfe("MH_Dfname=%s\n",MH_Dfname);
1.1       takayama   58:   }
                     59:   if (MH_RANK != mypower(2,MH_Mg)) {
1.13      takayama   60:     oxprintfe("rk.c MH_RANK=%d is not equal to 2^MH_Mg=%d. input_data_file is broken.\n",MH_RANK,mypower(2,MH_Mg));
1.7       takayama   61:     mh_usage();
                     62:     mh_exit(-1);
1.1       takayama   63:   }
                     64:   if (MH_Dfname != NULL) {
1.2       takayama   65:     if (!MH_RawName) {
                     66:       sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"-m%d-n%lf-x0%lf-b",MH_Mg,*MH_Ng,MH_X0g);
1.1       takayama   67:       for (i=0; i<MH_Mg; i++) {
1.7       takayama   68:         sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"%lf,",MH_Beta[i]);
1.1       takayama   69:       }
                     70:       sprintf(&(MH_Dfname[strlen(MH_Dfname)]),".txt");
                     71:     }
1.4       takayama   72:     Df = mh_fopen(MH_Dfname,"w",MH_byFile);
1.1       takayama   73:     if (Df == NULL) {
1.13      takayama   74:       oxprintfe("Error to open the file %s\n",MH_Dfname);
1.1       takayama   75:       return(-1);
                     76:     }
                     77:   }
                     78:   if (MH_Gfname != NULL) {
1.4       takayama   79:     Gf = mh_fopen(MH_Gfname,"w",MH_byFile);
1.1       takayama   80:     if (Gf == NULL) {
1.13      takayama   81:       oxprintfe("Error to open the file %s\n",MH_Gfname);
1.1       takayama   82:       return(-1);
                     83:     }
                     84:   }else return(1);
1.4       takayama   85:   if (MH_byFile) {
1.7       takayama   86:     sprintf(fname,"%s-gp.txt",MH_Gfname);
                     87:     fp = fopen(fname,"w");
                     88:     fprintf(fp,"set xrange [0:20]\n");
                     89:     fprintf(fp,"set yrange [0:1.2]\n");
                     90:     fprintf(fp,"plot \"%s\" title \"by hg\" with lines\n",MH_Gfname);
                     91:     fclose(fp);
1.4       takayama   92:   }
1.1       takayama   93:   return(0);
                     94: }
                     95: /*
1.7       takayama   96:   Runge-Kutta
                     97:   y_{n+1} = y_n + 1/6 k_1 + 1/3 k_2 + 1/3 k_3 + 1/6 k_4
                     98:   k_1 = h f(x_n, y_n)
                     99:   k_2 = h f(x_n + 1/2 h, y_n + 1/2 k_1)
                    100:   k_3 = h f(x_n + 1/2 h, y_n + 1/2 k_2)
                    101:   k_4 = h f(x_n + h, y_n + k_3)
1.1       takayama  102: */
                    103:
                    104: /* ベクトル値関数 rk is defined in tmp-code-?.c,
1.7       takayama  105:    f(x, y) (x : scalar, y : vector) */
1.1       takayama  106:
1.2       takayama  107: static void show_v(double x,double *v, int n)
1.1       takayama  108: {
                    109:   int i;
                    110:   static int counter=0;
1.2       takayama  111:   extern int MH_Dp;
1.4       takayama  112:   char swork[MH_SSIZE];
1.1       takayama  113:
1.16    ! takayama  114:   if (MH_Dp <= 0) return;
1.2       takayama  115:   if ((counter % MH_Dp) != 0) { counter++; return;} else counter=1;
1.4       takayama  116:   sprintf(swork,"%lf\n",x); mh_fputs(swork,Df);
                    117:   for (i = 0; i < n; i++) {sprintf(swork," %le\n", v[i]); mh_fputs(swork,Df);}
1.1       takayama  118: }
                    119:
1.4       takayama  120: struct MH_RESULT mh_rkmain(double x0,double y0[],double xn)
1.1       takayama  121: {
                    122:   static int initialized=0;
                    123:   int i;
                    124:   double h;
1.3       takayama  125:   double x;
1.4       takayama  126:   char swork[MH_SSIZE];
                    127:   struct MH_RESULT result;
1.3       takayama  128:   extern int MH_deallocate;
1.1       takayama  129:   /*
1.7       takayama  130:     double y[MH_RANK];
                    131:     double k1[MH_RANK], k2[MH_RANK], k3[MH_RANK], k4[MH_RANK];
                    132:     double temp[MH_RANK];
                    133:     double ty[MH_RANK];
1.1       takayama  134:   */
                    135:   static double *y,*k1,*k2,*k3,*k4,*temp,*ty;
1.14      takayama  136:
1.3       takayama  137:   if (MH_deallocate && initialized) {
1.7       takayama  138:     if (y) mh_free(y);
                    139:     if (k1) mh_free(k1);
                    140:     if (k2) mh_free(k2);
                    141:     if (k3) mh_free(k3);
                    142:     if (k4) mh_free(k4);
                    143:     if (temp) mh_free(temp);
                    144:     if (ty) mh_free(ty);
1.3       takayama  145:     y = k1 = k2 = k3 = k4 = temp = ty = NULL;
1.7       takayama  146:     initialized=0;
                    147:     return(result);
1.3       takayama  148:   }
1.1       takayama  149:   if (!initialized) {
1.2       takayama  150:     y = (double *)mh_malloc(sizeof(double)*MH_RANK);
                    151:     k1 = (double *)mh_malloc(sizeof(double)*MH_RANK);
                    152:     k2 = (double *)mh_malloc(sizeof(double)*MH_RANK);
                    153:     k3 = (double *)mh_malloc(sizeof(double)*MH_RANK);
                    154:     k4 = (double *)mh_malloc(sizeof(double)*MH_RANK);
                    155:     temp = (double *)mh_malloc(sizeof(double)*MH_RANK);
                    156:     ty = (double *)mh_malloc(sizeof(double)*MH_RANK);
1.7       takayama  157:     initialized=1;
1.1       takayama  158:   }
1.2       takayama  159:   h = MH_Hg;
1.1       takayama  160:   for (i = 0; i < MH_RANK; i++)
                    161:     y[i] = y0[i];
1.11      takayama  162: if (MH_strategy == 0) {
1.1       takayama  163:   for (x = x0; (h>0?(x<xn):(x>xn)); x += h) {
1.7       takayama  164:     if (Df) show_v(x,y, MH_RANK);
                    165:     if (Gf) {
                    166:       sprintf(swork,"%lf %le\n",x,y[0]);
                    167:       mh_fputs(swork,Gf);
                    168:     }
                    169:     /* Output 95% point */
                    170:     if (MH_P95) {
                    171:       if ((MH_P95==1) && (y[0] >= 0.95)) {
1.13      takayama  172:         oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
1.7       takayama  173:         MH_P95=2;
                    174:       }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {
1.13      takayama  175:         oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
1.7       takayama  176:         MH_P95=0;
                    177:       }
                    178:     }
1.15      takayama  179:     (*mh_rf)(x, y, MH_RANK, temp, MH_RANK);
1.1       takayama  180:     for (i = 0; i < MH_RANK; i++)
                    181:       k1[i] = h * temp[i];
                    182:
                    183:     for (i = 0; i < MH_RANK; i++)
                    184:       ty[i] = y[i] + 0.5 * k1[i];
1.15      takayama  185:     (*mh_rf)(x + 0.5 * h, ty, MH_RANK, temp, MH_RANK);
1.1       takayama  186:     for (i = 0; i < MH_RANK; i++)
                    187:       k2[i] = h * temp[i];
                    188:
                    189:     for (i = 0; i < MH_RANK; i++)
                    190:       ty[i] = y[i] + 0.5 * k2[i];
1.15      takayama  191:     (*mh_rf)(x + 0.5 * h, ty, MH_RANK, temp, MH_RANK);
1.1       takayama  192:     for (i = 0; i < MH_RANK; i++)
                    193:       k3[i] = h * temp[i];
                    194:
                    195:     for (i = 0; i < MH_RANK; i++)
                    196:       ty[i] = y[i] + k3[i];
1.15      takayama  197:     (*mh_rf)(x + h, ty, MH_RANK, temp, MH_RANK);
1.1       takayama  198:     for (i = 0; i < MH_RANK; i++)
                    199:       k4[i] = h * temp[i];
                    200:
                    201:     for (i = 0; i < MH_RANK; i++)
                    202:       y[i] = y[i] + 1.0/6.0 * k1[i] + 1.0/3.0 * k2[i] + 1.0/3.0 * k3[i] + 1.0/6.0 * k4[i];
                    203:   }
1.11      takayama  204: }else
1.9       takayama  205:   {
                    206:     extern int MH_Dp;
                    207:     double dh;
                    208:     double mh_dp_orig;
                    209:     double x1;
                    210:     const gsl_odeiv_step_type *T = gsl_odeiv_step_rkf45;
                    211:     gsl_odeiv_step *s = gsl_odeiv_step_alloc(T, MH_RANK);
1.11      takayama  212:     gsl_odeiv_control *c = gsl_odeiv_control_y_new(MH_abserr, MH_relerr);
1.9       takayama  213:     /* We should use the relative error.
                    214:        hgm.cwishart(m=5,n=20,beta=c(1,2,3,4,5),x=20,x0=2,approxdeg=20);
                    215:        0.977
                    216:        hgm.cwishart(m=8,n=20,beta=c(1,2,3,4,5,5.6,5.7,5.8),x=30,x0=1,approxdeg=20);
1.10      takayama  217:        0.962 (non-gsl) or 0.969 (gsl, abs error 0.0, relative error 1e-10),
                    218:        NaN abs=0, rel=1e-5 or abs=1e-18, rel=1e-6
                    219:        ./mh causes NaN when rel error is 1e-10 on old 32bit machines.
1.9       takayama  220:      */
                    221:     gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc(MH_RANK);
                    222:     gsl_odeiv_system sys = {mh_rf_for_gsl, NULL, 0, NULL};
                    223:     sys.dimension = MH_RANK;
1.13      takayama  224:     /* oxprintf("MH_RANK=%d\n",MH_RANK); */
                    225:     if (x0 >= xn) {oxprintfe("Error: x0 < x must hold.\n"); mh_exit(-30);}
1.9       takayama  226:     x = x0;
                    227:     if (MH_Dp > 0) dh = MH_Dp*h; else dh=xn-x0;
                    228:     mh_dp_orig = MH_Dp; MH_Dp=1;
                    229:     while (x < xn)  {
                    230:       if (Df) show_v(x,y, MH_RANK);
                    231:       if (Gf) {
                    232:         sprintf(swork,"%lf %le\n",x,y[0]);
                    233:         mh_fputs(swork,Gf);
                    234:       }
                    235:       /* Output 95% point */
                    236:       if (MH_P95) {
                    237:         if ((MH_P95==1) && (y[0] >= 0.95)) {
1.13      takayama  238:           oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
1.9       takayama  239:           MH_P95=2;
                    240:         }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {
1.13      takayama  241:           oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
1.9       takayama  242:           MH_P95=0;
                    243:         }
                    244:       }
                    245:       x1 = x+dh;
                    246:       while ((x < x1) && (x < xn)) {
                    247:         int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &x, x1, &h, y);
1.10      takayama  248:         if (status != GSL_SUCCESS) {
1.13      takayama  249:          oxprintfe("gsl_odeiv_evolve_apply failed.\n");
1.10      takayama  250:          break;
                    251:        }
1.9       takayama  252:       }
                    253:     }
                    254:     gsl_odeiv_evolve_free(e);
                    255:     gsl_odeiv_control_free(c);
                    256:     gsl_odeiv_step_free(s);
                    257:     MH_Dp=mh_dp_orig;
                    258:   }
1.11      takayama  259:
1.13      takayama  260:   if (MH_Verbose) oxprintf("x=%lf, y[0]=%lg\n",x,y[0]);
1.4       takayama  261:   result.x = x;
                    262:   result.rank = MH_RANK;
                    263:   result.y = (double *)mh_malloc(sizeof(double)*MH_RANK); /* todo, how to free it */
                    264:   for (i=0; i<MH_RANK; i++) (result.y)[i] = y[i];
                    265:   result.size=2;
                    266:   result.sfpp = (struct SFILE **)mh_malloc(sizeof(struct SFILE *)*(result.size)); /* todo, free */
                    267:   (result.sfpp)[0] = Df;
                    268:   (result.sfpp)[1] = Gf;
                    269:   return result;
1.1       takayama  270: }
1.11      takayama  271:
                    272:
                    273: /*
                    274:  rk4  MH_strategy==0;
                    275:  a-rk4  MH_strategy==1;  adaptive
                    276:  a-rk4-m  MH_strategy==2;  adaptive and multiply
                    277: */
                    278: void mh_set_strategy(int s,double err[2]) {
                    279:   MH_strategy = s;
                    280:   if (err[0] >= 0.0) MH_abserr = err[0];
                    281:   if (err[1] >= 0.0) MH_relerr = err[1];
                    282: }