[BACK]Return to rk.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / hgm / mh / src

Diff for /OpenXM/src/hgm/mh/src/rk.c between version 1.7 and 1.13

version 1.7, 2013/03/07 03:00:43 version 1.13, 2015/03/24 05:59:43
Line 1 
Line 1 
 /*  /*
   License: LGPL    License: LGPL
   Ref: Copied from this11/misc-2011/A1/wishart/Prog    Ref: Copied from this11/misc-2011/A1/wishart/Prog
   $OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.6 2013/03/01 05:26:25 takayama Exp $    $OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.12 2014/03/20 09:37:16 takayama Exp $
 */  */
 #include <stdio.h>  #include <stdio.h>
 #include <stdlib.h>  #include <stdlib.h>
 #include <math.h>  #include <math.h>
 #include <string.h>  #include <string.h>
 #include "sfile.h"  #include "sfile.h"
   #include "mh.h"
   
   #include "t-gsl_errno.h"
   #include "t-gsl_odeiv.h"
 extern int MH_RANK;  extern int MH_RANK;
 extern int MH_M;  extern int MH_M;
   
Line 22  static struct SFILE *Gf =NULL;
Line 25  static struct SFILE *Gf =NULL;
 static struct SFILE *Df =NULL;  static struct SFILE *Df =NULL;
 int MH_P95=0;  int MH_P95=0;
 int MH_Verbose=0;  int MH_Verbose=0;
   int MH_strategy=STRATEGY_DEFAULT;  /* 0: rk4, 1:a-rk4 (adaptive) */
   double MH_abserr = 1e-10;
   double MH_relerr = MH_RELERR_DEFAULT ;
   /* D_i = MH_abserr + MH_relerr*|y_i|
      If observed error E > D*110/100, the stepsize will be decreased.
      cf. GSL, 26.3 Adaptive Step-size control.
   */
   
   int mh_rf_for_gsl(double t,const double y[],double f[],void *params) {
     mh_rf(t,(double *)y,MH_RANK,f,MH_RANK);
     return(GSL_SUCCESS);
   }
   
 static int mypower(int x,int n) {  static int mypower(int x,int n) {
   int a,i;    int a,i;
   a = 1;    a = 1;
Line 36  int mh_gopen_file() {
Line 52  int mh_gopen_file() {
   Gf=NULL;    Gf=NULL;
   Df=NULL;    Df=NULL;
   if (MH_Verbose) {    if (MH_Verbose) {
     fprintf(stderr,"MH_Gfname=%s\n",MH_Gfname);      oxprintfe("MH_Gfname=%s\n",MH_Gfname);
     fprintf(stderr,"MH_Dfname=%s\n",MH_Dfname);      oxprintfe("MH_Dfname=%s\n",MH_Dfname);
   }    }
   if (MH_RANK != mypower(2,MH_Mg)) {    if (MH_RANK != mypower(2,MH_Mg)) {
     fprintf(stderr,"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));      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));
     mh_usage();      mh_usage();
     mh_exit(-1);      mh_exit(-1);
   }    }
Line 54  int mh_gopen_file() {
Line 70  int mh_gopen_file() {
     }      }
     Df = mh_fopen(MH_Dfname,"w",MH_byFile);      Df = mh_fopen(MH_Dfname,"w",MH_byFile);
     if (Df == NULL) {      if (Df == NULL) {
       fprintf(stderr,"Error to open the file %s\n",MH_Dfname);        oxprintfe("Error to open the file %s\n",MH_Dfname);
       return(-1);        return(-1);
     }      }
   }    }
   if (MH_Gfname != NULL) {    if (MH_Gfname != NULL) {
     Gf = mh_fopen(MH_Gfname,"w",MH_byFile);      Gf = mh_fopen(MH_Gfname,"w",MH_byFile);
     if (Gf == NULL) {      if (Gf == NULL) {
       fprintf(stderr,"Error to open the file %s\n",MH_Gfname);        oxprintfe("Error to open the file %s\n",MH_Gfname);
       return(-1);        return(-1);
     }      }
   }else return(1);    }else return(1);
Line 140  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
Line 156  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
   h = MH_Hg;    h = MH_Hg;
   for (i = 0; i < MH_RANK; i++)    for (i = 0; i < MH_RANK; i++)
     y[i] = y0[i];      y[i] = y0[i];
   if (MH_strategy == 0) {
   for (x = x0; (h>0?(x<xn):(x>xn)); x += h) {    for (x = x0; (h>0?(x<xn):(x>xn)); x += h) {
     if (Df) show_v(x,y, MH_RANK);      if (Df) show_v(x,y, MH_RANK);
     if (Gf) {      if (Gf) {
Line 149  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
Line 166  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
     /* Output 95% point */      /* Output 95% point */
     if (MH_P95) {      if (MH_P95) {
       if ((MH_P95==1) && (y[0] >= 0.95)) {        if ((MH_P95==1) && (y[0] >= 0.95)) {
         printf("x=%le, y[0]=%lf\n",x,y[0]);          oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
         MH_P95=2;          MH_P95=2;
       }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {        }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {
         printf("x=%le, y[0]=%lf\n",x,y[0]);          oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
         MH_P95=0;          MH_P95=0;
       }        }
     }      }
Line 181  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
Line 198  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
     for (i = 0; i < MH_RANK; i++)      for (i = 0; i < MH_RANK; i++)
       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];        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];
   }    }
   if (MH_Verbose) printf("x=%lf, y[0]=%lg\n",x,y[0]);  }else
     {
       extern int MH_Dp;
       double dh;
       double mh_dp_orig;
       double x1;
       const gsl_odeiv_step_type *T = gsl_odeiv_step_rkf45;
       gsl_odeiv_step *s = gsl_odeiv_step_alloc(T, MH_RANK);
       gsl_odeiv_control *c = gsl_odeiv_control_y_new(MH_abserr, MH_relerr);
       /* We should use the relative error.
          hgm.cwishart(m=5,n=20,beta=c(1,2,3,4,5),x=20,x0=2,approxdeg=20);
          0.977
          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);
          0.962 (non-gsl) or 0.969 (gsl, abs error 0.0, relative error 1e-10),
          NaN abs=0, rel=1e-5 or abs=1e-18, rel=1e-6
          ./mh causes NaN when rel error is 1e-10 on old 32bit machines.
        */
       gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc(MH_RANK);
       gsl_odeiv_system sys = {mh_rf_for_gsl, NULL, 0, NULL};
       sys.dimension = MH_RANK;
       /* oxprintf("MH_RANK=%d\n",MH_RANK); */
       if (x0 >= xn) {oxprintfe("Error: x0 < x must hold.\n"); mh_exit(-30);}
       x = x0;
       if (MH_Dp > 0) dh = MH_Dp*h; else dh=xn-x0;
       mh_dp_orig = MH_Dp; MH_Dp=1;
       while (x < xn)  {
         if (Df) show_v(x,y, MH_RANK);
         if (Gf) {
           sprintf(swork,"%lf %le\n",x,y[0]);
           mh_fputs(swork,Gf);
         }
         /* Output 95% point */
         if (MH_P95) {
           if ((MH_P95==1) && (y[0] >= 0.95)) {
             oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
             MH_P95=2;
           }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {
             oxprintf("x=%le, y[0]=%lf\n",x,y[0]);
             MH_P95=0;
           }
         }
         x1 = x+dh;
         while ((x < x1) && (x < xn)) {
           int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &x, x1, &h, y);
           if (status != GSL_SUCCESS) {
             oxprintfe("gsl_odeiv_evolve_apply failed.\n");
             break;
           }
         }
       }
       gsl_odeiv_evolve_free(e);
       gsl_odeiv_control_free(c);
       gsl_odeiv_step_free(s);
       MH_Dp=mh_dp_orig;
     }
   
     if (MH_Verbose) oxprintf("x=%lf, y[0]=%lg\n",x,y[0]);
   result.x = x;    result.x = x;
   result.rank = MH_RANK;    result.rank = MH_RANK;
   result.y = (double *)mh_malloc(sizeof(double)*MH_RANK); /* todo, how to free it */    result.y = (double *)mh_malloc(sizeof(double)*MH_RANK); /* todo, how to free it */
Line 191  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
Line 264  struct MH_RESULT mh_rkmain(double x0,double y0[],doubl
   (result.sfpp)[0] = Df;    (result.sfpp)[0] = Df;
   (result.sfpp)[1] = Gf;    (result.sfpp)[1] = Gf;
   return result;    return result;
   }
   
   
   /*
    rk4  MH_strategy==0;
    a-rk4  MH_strategy==1;  adaptive
    a-rk4-m  MH_strategy==2;  adaptive and multiply
   */
   void mh_set_strategy(int s,double err[2]) {
     MH_strategy = s;
     if (err[0] >= 0.0) MH_abserr = err[0];
     if (err[1] >= 0.0) MH_relerr = err[1];
 }  }
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.13
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>