[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.3 and 1.4

version 1.3, 2013/02/20 01:06:38 version 1.4, 2013/02/20 05:20:49
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.2 2013/02/19 08:03:14 takayama Exp $  $OpenXM: OpenXM/src/hgm/mh/src/rk.c,v 1.3 2013/02/20 01:06:38 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"
   
 extern int MH_RANK;  extern int MH_RANK;
 extern int MH_M;  extern int MH_M;
Line 17  extern int MH_Mg;
Line 18  extern int MH_Mg;
 extern double *MH_Beta; extern double *MH_Ng; extern double MH_X0g;  extern double *MH_Beta; extern double *MH_Ng; extern double MH_X0g;
 extern int MH_RawName;  extern int MH_RawName;
 double MH_Hg = 0.001;  double MH_Hg = 0.001;
 static FILE *Gf=NULL;  static struct SFILE *Gf =NULL;
 static FILE *Df =NULL;  static struct SFILE *Df =NULL;
 int MH_P95=0;  int MH_P95=0;
 int MH_Verbose=0;  int MH_Verbose=0;
 static mypower(int x,int n) {  static mypower(int x,int n) {
Line 31  mh_gopen_file() {
Line 32  mh_gopen_file() {
   FILE *fp;    FILE *fp;
   char fname[1024];    char fname[1024];
   int i;    int i;
     extern int MH_byFile;
   Gf=NULL;    Gf=NULL;
   Df=NULL;    Df=NULL;
   if (MH_Verbose) {    if (MH_Verbose) {
Line 46  mh_gopen_file() {
Line 48  mh_gopen_file() {
     if (!MH_RawName) {      if (!MH_RawName) {
       sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"-m%d-n%lf-x0%lf-b",MH_Mg,*MH_Ng,MH_X0g);        sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"-m%d-n%lf-x0%lf-b",MH_Mg,*MH_Ng,MH_X0g);
       for (i=0; i<MH_Mg; i++) {        for (i=0; i<MH_Mg; i++) {
         sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"%lf,",MH_Beta[i]);                  sprintf(&(MH_Dfname[strlen(MH_Dfname)]),"%lf,",MH_Beta[i]);
       }        }
       sprintf(&(MH_Dfname[strlen(MH_Dfname)]),".txt");        sprintf(&(MH_Dfname[strlen(MH_Dfname)]),".txt");
     }      }
     Df = fopen(MH_Dfname,"w");      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);        fprintf(stderr,"Error to open the file %s\n",MH_Dfname);
       return(-1);        return(-1);
     }      }
   }    }
   if (MH_Gfname != NULL) {    if (MH_Gfname != NULL) {
     Gf = fopen(MH_Gfname,"w");      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);        fprintf(stderr,"Error to open the file %s\n",MH_Gfname);
       return(-1);        return(-1);
     }      }
   }else return(1);    }else return(1);
   sprintf(fname,"%s-gp.txt",MH_Gfname);    if (MH_byFile) {
   fp = fopen(fname,"w");          sprintf(fname,"%s-gp.txt",MH_Gfname);
   fprintf(fp,"set xrange [0:20]\n");          fp = fopen(fname,"w");
   fprintf(fp,"set yrange [0:1.2]\n");          fprintf(fp,"set xrange [0:20]\n");
   fprintf(fp,"plot \"%s\" title \"by hg\" with lines\n",MH_Gfname);          fprintf(fp,"set yrange [0:1.2]\n");
   fclose(fp);          fprintf(fp,"plot \"%s\" title \"by hg\" with lines\n",MH_Gfname);
           fclose(fp);
     }
   return(0);    return(0);
 }  }
 /*  /*
Line 88  static void show_v(double x,double *v, int n)
Line 92  static void show_v(double x,double *v, int n)
   int i;    int i;
   static int counter=0;    static int counter=0;
   extern int MH_Dp;    extern int MH_Dp;
     char swork[MH_SSIZE];
   
   if ((counter % MH_Dp) != 0) { counter++; return;} else counter=1;    if ((counter % MH_Dp) != 0) { counter++; return;} else counter=1;
   fprintf(Df,"%lf\n",x);    sprintf(swork,"%lf\n",x); mh_fputs(swork,Df);
   for (i = 0; i < n; i++) fprintf(Df," %le\n", v[i]);    for (i = 0; i < n; i++) {sprintf(swork," %le\n", v[i]); mh_fputs(swork,Df);}
 }  }
   
 double mh_rkmain(double x0,double y0[],double xn)  struct MH_RESULT mh_rkmain(double x0,double y0[],double xn)
 {  {
   static int initialized=0;    static int initialized=0;
   int i;    int i;
   double h;    double h;
   double x;    double x;
     char swork[MH_SSIZE];
     struct MH_RESULT result;
   extern int MH_deallocate;    extern int MH_deallocate;
   /*    /*
   double y[MH_RANK];    double y[MH_RANK];
Line 118  double mh_rkmain(double x0,double y0[],double xn)
Line 125  double mh_rkmain(double x0,double y0[],double xn)
         if (ty) mh_free(ty);          if (ty) mh_free(ty);
     y = k1 = k2 = k3 = k4 = temp = ty = NULL;      y = k1 = k2 = k3 = k4 = temp = ty = NULL;
         initialized=0;          initialized=0;
         return(0.0);          return(result);
   }    }
   if (!initialized) {    if (!initialized) {
     y = (double *)mh_malloc(sizeof(double)*MH_RANK);      y = (double *)mh_malloc(sizeof(double)*MH_RANK);
Line 135  double mh_rkmain(double x0,double y0[],double xn)
Line 142  double mh_rkmain(double x0,double y0[],double xn)
     y[i] = y0[i];      y[i] = y0[i];
   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) fprintf(Gf,"%lf %le\n",x,y[0]);          if (Gf) {
             sprintf(swork,"%lf %le\n",x,y[0]);
             mh_fputs(swork,Gf);
           }
         /* 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)) {
Line 172  double mh_rkmain(double x0,double y0[],double xn)
Line 182  double mh_rkmain(double x0,double y0[],double xn)
       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];
   }    }
   printf("x=%lf, y[0]=%lg\n",x,y[0]);    printf("x=%lf, y[0]=%lg\n",x,y[0]);
   if (Df) fclose(Df);    result.x = x;
   if (Gf) fclose(Gf);    result.rank = MH_RANK;
   return x;    result.y = (double *)mh_malloc(sizeof(double)*MH_RANK); /* todo, how to free it */
     for (i=0; i<MH_RANK; i++) (result.y)[i] = y[i];
     result.size=2;
     result.sfpp = (struct SFILE **)mh_malloc(sizeof(struct SFILE *)*(result.size)); /* todo, free */
     (result.sfpp)[0] = Df;
     (result.sfpp)[1] = Gf;
     return result;
 }  }
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>