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

Diff for /OpenXM/src/hgm/mh/src/code-n.c between version 1.3 and 1.6

version 1.3, 2013/02/20 01:06:37 version 1.6, 2013/03/07 03:00:43
Line 1 
Line 1 
 /*  /*
 $OpenXM: OpenXM/src/hgm/mh/src/code-n.c,v 1.2 2013/02/19 08:03:14 takayama Exp $  $OpenXM: OpenXM/src/hgm/mh/src/code-n.c,v 1.5 2013/03/05 05:26:07 takayama Exp $
 License: LGPL  License: LGPL
 Ref: Copied from this11/misc-2011/A1/wishart/Prog  Ref: Copied from this11/misc-2011/A1/wishart/Prog
 cf. @s/2011/12/01-my-note-code-n.pdf  cf. @s/2011/12/01-my-note-code-n.pdf
Line 7  cf. @s/2011/12/01-my-note-code-n.pdf
Line 7  cf. @s/2011/12/01-my-note-code-n.pdf
 #include <stdio.h>  #include <stdio.h>
 #include <stdlib.h>  #include <stdlib.h>
 #include <math.h>  #include <math.h>
   #include "sfile.h"
   
 static void error_code(char *s);  static void error_code(char *s);
 static void showf(char *s,double *v,int m);  static void showf(char *s,double *v,int m);
Line 60  void mh_rf(double x, double *f, int rank_not_used, dou
Line 61  void mh_rf(double x, double *f, int rank_not_used, dou
   /* double f2[MH_M][MH_RANK];*/    /* double f2[MH_M][MH_RANK];*/
   static double *f2=NULL;    static double *f2=NULL;
   
     mh_check_intr(100);
   if (MH_deallocate && initialized) {    if (MH_deallocate && initialized) {
         if (b) mh_free(b);          if (b) mh_free(b);
         if (bitSize) mh_free(bitSize);          if (bitSize) mh_free(bitSize);
         if (invbitSize) mh_free(invbitSize);          if (invbitSize) mh_free(invbitSize);
         if (f2) mh_free(f2);          if (f2) mh_free(f2);
         if (yik) mh_free(yik);          if (yik) mh_free(yik);
         if (yik2) mh_free(yik2);          if (yik2) mh_free(yik2);
         if (y) mh_free(y);          if (y) mh_free(y);
     b = f2 = yik = yik2 = y = NULL;      b = f2 = yik = yik2 = y = NULL;
     bitSize = invbitSize = NULL;      bitSize = invbitSize = NULL;
         initialized = 0; return;          initialized = 0; return;
   }    }
   if (!initialized) {    if (!initialized) {
     b = (double *)mh_malloc(sizeof(double)*MH_M);      b = (double *)mh_malloc(sizeof(double)*MH_M);
Line 96  void mh_rf(double x, double *f, int rank_not_used, dou
Line 98  void mh_rf(double x, double *f, int rank_not_used, dou
     for (i=0; i<MH_RANK; i++) bitSize[i] = bitcount(i);      for (i=0; i<MH_RANK; i++) bitSize[i] = bitcount(i);
     for (i=0, k=0; i<=MH_M; i++) {      for (i=0, k=0; i<=MH_M; i++) {
       for (j=0; j<MH_RANK; j++) {        for (j=0; j<MH_RANK; j++) {
         if (bitSize[j] == i) {invbitSize[k] = j; k++; }          if (bitSize[j] == i) {invbitSize[k] = j; k++; }
       }        }
     }      }
     initialized = 1;      initialized = 1;
Line 107  void mh_rf(double x, double *f, int rank_not_used, dou
Line 109  void mh_rf(double x, double *f, int rank_not_used, dou
   for (i=0; i<MH_M; i++) {    for (i=0; i<MH_M; i++) {
     for (k=0; k<MH_M; k++) {      for (k=0; k<MH_M; k++) {
       if (i != k) {        if (i != k) {
         yik[idxM(i,k)] = y[k]/(y[i]-y[k]);          yik[idxM(i,k)] = y[k]/(y[i]-y[k]);
         yik2[idxM(i,k)] = y[i]/((y[i]-y[k])*(y[i]-y[k]));          yik2[idxM(i,k)] = y[i]/((y[i]-y[k])*(y[i]-y[k]));
       }        }
     }      }
   }    }
Line 122  void mh_rf(double x, double *f, int rank_not_used, dou
Line 124  void mh_rf(double x, double *f, int rank_not_used, dou
     f2[idxRank(i,0)] = -((c-y[i])*f[join(i,0)]-a*f[0]);      f2[idxRank(i,0)] = -((c-y[i])*f[join(i,0)]-a*f[0]);
     for (k=0; k<MH_M; k++) {      for (k=0; k<MH_M; k++) {
       if (i!=k) {        if (i!=k) {
         f2[idxRank(i,0)] += -0.5*yik[idxM(i,k)]*(f[join(i,0)]-f[join(k,0)]);          f2[idxRank(i,0)] += -0.5*yik[idxM(i,k)]*(f[join(i,0)]-f[join(k,0)]);
       }        }
     }      }
     f2[idxRank(i,0)] = f2[idxRank(i,0)]/y[i];      f2[idxRank(i,0)] = f2[idxRank(i,0)]/y[i];
Line 134  void mh_rf(double x, double *f, int rank_not_used, dou
Line 136  void mh_rf(double x, double *f, int rank_not_used, dou
       ii = join(i,jj);        ii = join(i,jj);
       rijj = -( (c-y[i])*f[ii]-a*f[jj] );        rijj = -( (c-y[i])*f[ii]-a*f[jj] );
       for (k=0; k<MH_M; k++) {        for (k=0; k<MH_M; k++) {
         if (!member(k,ii)) rijj += -(0.5)*yik[idxM(i,k)]*(f[ii]-f[join(k,jj)]);          if (!member(k,ii)) rijj += -(0.5)*yik[idxM(i,k)]*(f[ii]-f[join(k,jj)]);
         if (member(k,jj)) {          if (member(k,jj)) {
           rijj += -(0.5)*yik[idxM(i,k)]*f[ii];            rijj += -(0.5)*yik[idxM(i,k)]*f[ii];
           rijj += -(0.5)*yik2[idxM(i,k)]*(f[join(i,delete(k,jj))]-f[jj]);            rijj += -(0.5)*yik2[idxM(i,k)]*(f[join(i,delete(k,jj))]-f[jj]);
         }          }
       }        }
       f2[idxRank(i,jj)] = rijj;        f2[idxRank(i,jj)] = rijj;
       for (k = 0; k<MH_M; k++) {        for (k = 0; k<MH_M; k++) {
         if (member(k,jj)) {          if (member(k,jj)) {
           f2[idxRank(i,jj)] += (0.5)*yik[idxM(i,k)]*f2[idxRank(k,delete(k,jj))];            f2[idxRank(i,jj)] += (0.5)*yik[idxM(i,k)]*f2[idxRank(k,delete(k,jj))];
         }          }
       }        }
       f2[idxRank(i,jj)] = f2[idxRank(i,jj)]/y[i];        f2[idxRank(i,jj)] = f2[idxRank(i,jj)]/y[i];
     }      }
Line 155  void mh_rf(double x, double *f, int rank_not_used, dou
Line 157  void mh_rf(double x, double *f, int rank_not_used, dou
     val[jj] = 0;      val[jj] = 0;
     for (i=0; i<MH_M; i++) {      for (i=0; i<MH_M; i++) {
       if (member(i,jj)) {        if (member(i,jj)) {
         val[jj] += b[i]*f2[idxRank(i,delete(i,jj))];          val[jj] += b[i]*f2[idxRank(i,delete(i,jj))];
       }else{        }else{
         val[jj] += b[i]*f[join(i,jj)];          val[jj] += b[i]*f[join(i,jj)];
       }        }
     }      }
   }    }
Line 203  showf2(char *s,double *v,int m,int n) {
Line 205  showf2(char *s,double *v,int m,int n) {
   }    }
   printf("\n");    printf("\n");
   printf("member(0,0)=%d, member(0,1)=%d, member(0,2)=%d,member(0,3)=%d\n",    printf("member(0,0)=%d, member(0,1)=%d, member(0,2)=%d,member(0,3)=%d\n",
          member(0,0),member(0,1),member(0,2),member(0,3));           member(0,0),member(0,1),member(0,2),member(0,3));
 }  }
   
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>