Annotation of OpenXM/src/hgm/mh/src/rk.c, Revision 1.9
1.1 takayama 1: /*
1.7 takayama 2: License: LGPL
3: Ref: Copied from this11/misc-2011/A1/wishart/Prog
1.9 ! takayama 4: $OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.8 2013/03/07 05:23:31 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;
17:
18: char *MH_Gfname;
19: char *MH_Dfname;
20: extern int MH_Mg;
1.2 takayama 21: extern double *MH_Beta; extern double *MH_Ng; extern double MH_X0g;
22: extern int MH_RawName;
23: double MH_Hg = 0.001;
1.4 takayama 24: static struct SFILE *Gf =NULL;
25: static struct SFILE *Df =NULL;
1.2 takayama 26: int MH_P95=0;
27: int MH_Verbose=0;
1.9 ! takayama 28:
! 29: int mh_rf_for_gsl(double t,const double y[],double f[],void *params) {
! 30: mh_rf(t,(double *)y,MH_RANK,f,MH_RANK);
! 31: return(GSL_SUCCESS);
! 32: }
! 33:
1.6 takayama 34: static int mypower(int x,int n) {
1.1 takayama 35: int a,i;
36: a = 1;
37: for (i=0; i<n; i++) a = a*x;
38: return(a);
39: }
1.6 takayama 40: int mh_gopen_file() {
1.1 takayama 41: FILE *fp;
42: char fname[1024];
43: int i;
1.4 takayama 44: extern int MH_byFile;
1.1 takayama 45: Gf=NULL;
46: Df=NULL;
1.2 takayama 47: if (MH_Verbose) {
1.1 takayama 48: fprintf(stderr,"MH_Gfname=%s\n",MH_Gfname);
49: fprintf(stderr,"MH_Dfname=%s\n",MH_Dfname);
50: }
51: if (MH_RANK != mypower(2,MH_Mg)) {
1.7 takayama 52: 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));
53: mh_usage();
54: mh_exit(-1);
1.1 takayama 55: }
56: if (MH_Dfname != NULL) {
1.2 takayama 57: if (!MH_RawName) {
58: sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"-m%d-n%lf-x0%lf-b",MH_Mg,*MH_Ng,MH_X0g);
1.1 takayama 59: for (i=0; i<MH_Mg; i++) {
1.7 takayama 60: sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"%lf,",MH_Beta[i]);
1.1 takayama 61: }
62: sprintf(&(MH_Dfname[strlen(MH_Dfname)]),".txt");
63: }
1.4 takayama 64: Df = mh_fopen(MH_Dfname,"w",MH_byFile);
1.1 takayama 65: if (Df == NULL) {
66: fprintf(stderr,"Error to open the file %s\n",MH_Dfname);
67: return(-1);
68: }
69: }
70: if (MH_Gfname != NULL) {
1.4 takayama 71: Gf = mh_fopen(MH_Gfname,"w",MH_byFile);
1.1 takayama 72: if (Gf == NULL) {
73: fprintf(stderr,"Error to open the file %s\n",MH_Gfname);
74: return(-1);
75: }
76: }else return(1);
1.4 takayama 77: if (MH_byFile) {
1.7 takayama 78: sprintf(fname,"%s-gp.txt",MH_Gfname);
79: fp = fopen(fname,"w");
80: fprintf(fp,"set xrange [0:20]\n");
81: fprintf(fp,"set yrange [0:1.2]\n");
82: fprintf(fp,"plot \"%s\" title \"by hg\" with lines\n",MH_Gfname);
83: fclose(fp);
1.4 takayama 84: }
1.1 takayama 85: return(0);
86: }
87: /*
1.7 takayama 88: Runge-Kutta
89: y_{n+1} = y_n + 1/6 k_1 + 1/3 k_2 + 1/3 k_3 + 1/6 k_4
90: k_1 = h f(x_n, y_n)
91: k_2 = h f(x_n + 1/2 h, y_n + 1/2 k_1)
92: k_3 = h f(x_n + 1/2 h, y_n + 1/2 k_2)
93: k_4 = h f(x_n + h, y_n + k_3)
1.1 takayama 94: */
95:
96: /* ベクトル値関数 rk is defined in tmp-code-?.c,
1.7 takayama 97: f(x, y) (x : scalar, y : vector) */
1.1 takayama 98:
1.2 takayama 99: static void show_v(double x,double *v, int n)
1.1 takayama 100: {
101: int i;
102: static int counter=0;
1.2 takayama 103: extern int MH_Dp;
1.4 takayama 104: char swork[MH_SSIZE];
1.1 takayama 105:
1.2 takayama 106: if ((counter % MH_Dp) != 0) { counter++; return;} else counter=1;
1.4 takayama 107: sprintf(swork,"%lf\n",x); mh_fputs(swork,Df);
108: for (i = 0; i < n; i++) {sprintf(swork," %le\n", v[i]); mh_fputs(swork,Df);}
1.1 takayama 109: }
110:
1.4 takayama 111: struct MH_RESULT mh_rkmain(double x0,double y0[],double xn)
1.1 takayama 112: {
113: static int initialized=0;
114: int i;
115: double h;
1.3 takayama 116: double x;
1.4 takayama 117: char swork[MH_SSIZE];
118: struct MH_RESULT result;
1.3 takayama 119: extern int MH_deallocate;
1.1 takayama 120: /*
1.7 takayama 121: double y[MH_RANK];
122: double k1[MH_RANK], k2[MH_RANK], k3[MH_RANK], k4[MH_RANK];
123: double temp[MH_RANK];
124: double ty[MH_RANK];
1.1 takayama 125: */
126: static double *y,*k1,*k2,*k3,*k4,*temp,*ty;
1.3 takayama 127: if (MH_deallocate && initialized) {
1.7 takayama 128: if (y) mh_free(y);
129: if (k1) mh_free(k1);
130: if (k2) mh_free(k2);
131: if (k3) mh_free(k3);
132: if (k4) mh_free(k4);
133: if (temp) mh_free(temp);
134: if (ty) mh_free(ty);
1.3 takayama 135: y = k1 = k2 = k3 = k4 = temp = ty = NULL;
1.7 takayama 136: initialized=0;
137: return(result);
1.3 takayama 138: }
1.1 takayama 139: if (!initialized) {
1.2 takayama 140: y = (double *)mh_malloc(sizeof(double)*MH_RANK);
141: k1 = (double *)mh_malloc(sizeof(double)*MH_RANK);
142: k2 = (double *)mh_malloc(sizeof(double)*MH_RANK);
143: k3 = (double *)mh_malloc(sizeof(double)*MH_RANK);
144: k4 = (double *)mh_malloc(sizeof(double)*MH_RANK);
145: temp = (double *)mh_malloc(sizeof(double)*MH_RANK);
146: ty = (double *)mh_malloc(sizeof(double)*MH_RANK);
1.7 takayama 147: initialized=1;
1.1 takayama 148: }
1.2 takayama 149: h = MH_Hg;
1.1 takayama 150: for (i = 0; i < MH_RANK; i++)
151: y[i] = y0[i];
1.9 ! takayama 152: #ifdef NO_GSL
1.1 takayama 153: for (x = x0; (h>0?(x<xn):(x>xn)); x += h) {
1.7 takayama 154: if (Df) show_v(x,y, MH_RANK);
155: if (Gf) {
156: sprintf(swork,"%lf %le\n",x,y[0]);
157: mh_fputs(swork,Gf);
158: }
159: /* Output 95% point */
160: if (MH_P95) {
161: if ((MH_P95==1) && (y[0] >= 0.95)) {
162: printf("x=%le, y[0]=%lf\n",x,y[0]);
163: MH_P95=2;
164: }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {
165: printf("x=%le, y[0]=%lf\n",x,y[0]);
166: MH_P95=0;
167: }
168: }
1.2 takayama 169: mh_rf(x, y, MH_RANK, temp, MH_RANK);
1.1 takayama 170: for (i = 0; i < MH_RANK; i++)
171: k1[i] = h * temp[i];
172:
173: for (i = 0; i < MH_RANK; i++)
174: ty[i] = y[i] + 0.5 * k1[i];
1.2 takayama 175: mh_rf(x + 0.5 * h, ty, MH_RANK, temp, MH_RANK);
1.1 takayama 176: for (i = 0; i < MH_RANK; i++)
177: k2[i] = h * temp[i];
178:
179: for (i = 0; i < MH_RANK; i++)
180: ty[i] = y[i] + 0.5 * k2[i];
1.2 takayama 181: mh_rf(x + 0.5 * h, ty, MH_RANK, temp, MH_RANK);
1.1 takayama 182: for (i = 0; i < MH_RANK; i++)
183: k3[i] = h * temp[i];
184:
185: for (i = 0; i < MH_RANK; i++)
186: ty[i] = y[i] + k3[i];
1.2 takayama 187: mh_rf(x + h, ty, MH_RANK, temp, MH_RANK);
1.1 takayama 188: for (i = 0; i < MH_RANK; i++)
189: k4[i] = h * temp[i];
190:
191: for (i = 0; i < MH_RANK; i++)
192: 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];
193: }
1.9 ! takayama 194: #else
! 195: {
! 196: extern int MH_Dp;
! 197: double dh;
! 198: double mh_dp_orig;
! 199: double x1;
! 200: const gsl_odeiv_step_type *T = gsl_odeiv_step_rkf45;
! 201: gsl_odeiv_step *s = gsl_odeiv_step_alloc(T, MH_RANK);
! 202: gsl_odeiv_control *c = gsl_odeiv_control_y_new(0.0, 1e-10);
! 203: /* We should use the relative error.
! 204: hgm.cwishart(m=5,n=20,beta=c(1,2,3,4,5),x=20,x0=2,approxdeg=20);
! 205: 0.977
! 206: 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);
! 207: 0.962 (non-gsl) or 0.969 (gsl)
! 208: */
! 209: gsl_odeiv_evolve *e = gsl_odeiv_evolve_alloc(MH_RANK);
! 210: gsl_odeiv_system sys = {mh_rf_for_gsl, NULL, 0, NULL};
! 211: sys.dimension = MH_RANK;
! 212: /* printf("MH_RANK=%d\n",MH_RANK); */
! 213: if (x0 >= xn) {fprintf(stderr,"Error: x0 < x must hold.\n"); mh_exit(-30);}
! 214: x = x0;
! 215: if (MH_Dp > 0) dh = MH_Dp*h; else dh=xn-x0;
! 216: mh_dp_orig = MH_Dp; MH_Dp=1;
! 217: while (x < xn) {
! 218: if (Df) show_v(x,y, MH_RANK);
! 219: if (Gf) {
! 220: sprintf(swork,"%lf %le\n",x,y[0]);
! 221: mh_fputs(swork,Gf);
! 222: }
! 223: /* Output 95% point */
! 224: if (MH_P95) {
! 225: if ((MH_P95==1) && (y[0] >= 0.95)) {
! 226: printf("x=%le, y[0]=%lf\n",x,y[0]);
! 227: MH_P95=2;
! 228: }else if ((MH_P95==2) && (y[0] >=0.9500000001)) {
! 229: printf("x=%le, y[0]=%lf\n",x,y[0]);
! 230: MH_P95=0;
! 231: }
! 232: }
! 233: x1 = x+dh;
! 234: while ((x < x1) && (x < xn)) {
! 235: int status = gsl_odeiv_evolve_apply(e, c, s, &sys, &x, x1, &h, y);
! 236: if (status != GSL_SUCCESS) break;
! 237: }
! 238: }
! 239: gsl_odeiv_evolve_free(e);
! 240: gsl_odeiv_control_free(c);
! 241: gsl_odeiv_step_free(s);
! 242: MH_Dp=mh_dp_orig;
! 243: }
! 244: #endif
1.5 takayama 245: if (MH_Verbose) printf("x=%lf, y[0]=%lg\n",x,y[0]);
1.4 takayama 246: result.x = x;
247: result.rank = MH_RANK;
248: result.y = (double *)mh_malloc(sizeof(double)*MH_RANK); /* todo, how to free it */
249: for (i=0; i<MH_RANK; i++) (result.y)[i] = y[i];
250: result.size=2;
251: result.sfpp = (struct SFILE **)mh_malloc(sizeof(struct SFILE *)*(result.size)); /* todo, free */
252: (result.sfpp)[0] = Df;
253: (result.sfpp)[1] = Gf;
254: return result;
1.1 takayama 255: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to rk.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / hgm / mh / src