Annotation of OpenXM/src/hgm/mh/src/code-n.c, Revision 1.1
1.1 ! takayama 1: /*
! 2: License: LGPL
! 3: cf. @s/2011/12/01-my-note-code-n.pdf
! 4: */
! 5: #include <stdio.h>
! 6: #include <stdlib.h>
! 7: #include <math.h>
! 8:
! 9: extern int MH_M;
! 10: extern int MH_RANK;
! 11:
! 12: int bitcount(int n) {int c; c=0; while (n) { ++c; n &= (n-1); } return(c);}
! 13: #define join(k,jj) ( (1 << k) | jj)
! 14: #define delete(k,jj) ( (~(1 << k)) & jj)
! 15: #define member(k,ii) ( (1 << k) & ii)
! 16: #define NaN 9.999e100
! 17:
! 18: #define idxM(i,j) ((i)*MH_M+j)
! 19: #define idxRank(i,j) ((i)*MH_RANK+j)
! 20:
! 21: /*
! 22: p(x) = exp(-bsum*x)*x^(MH_M*(c-a))
! 23: f is p(x)*pd{J}F
! 24: */
! 25: void rf(double x, double *f, int rank_not_used, double *val, int n_not_used)
! 26: { extern double *MH_Beta; /* beta = (1/2)*\simga^(-1) */
! 27: extern double *Ng; /* freedom */
! 28: extern int MH_Mg; /* number of variables, MH_Mg=MH_M */
! 29: static int initialized = 0;
! 30: /* static double b[MH_M]; */
! 31: static double *b;
! 32: static double bsum = 0;
! 33: static double a;
! 34: static double c;
! 35: static int m;
! 36: /*
! 37: static int bitSize[MH_RANK];
! 38: static int invbitSize[MH_RANK];
! 39: */
! 40: static int *bitSize;
! 41: static int *invbitSize;
! 42: int i,j,k,p;
! 43: int ii,jj; /* stands for I and J */
! 44: double rijj;
! 45: /*
! 46: double yik[MH_M][MH_M];
! 47: double yik2[MH_M][MH_M];
! 48: double y[MH_M];
! 49: */
! 50: static double *yik=NULL;
! 51: static double *yik2=NULL;
! 52: static double *y=NULL;
! 53: /* double f2[MH_M][MH_RANK];*/
! 54: static double *f2=NULL;
! 55:
! 56: if (!initialized) {
! 57: b = (double *)mymalloc(sizeof(double)*MH_M);
! 58: bitSize = (int *)mymalloc(sizeof(int)*MH_RANK);
! 59: invbitSize = (int *)mymalloc(sizeof(int)*MH_RANK);
! 60:
! 61: f2 = (double *)mymalloc(sizeof(double)*MH_M*MH_RANK);
! 62: yik = (double *)mymalloc(sizeof(double)*MH_M*MH_M);
! 63: yik2 = (double *)mymalloc(sizeof(double)*MH_M*MH_M);
! 64: y = (double *)mymalloc(sizeof(double)*MH_M);
! 65:
! 66: m = MH_Mg;
! 67: if (m != MH_M) error_code("MH_M != m. MH_M is given by -DMH_M...");
! 68: a = ((double)m+1.0)/2.0;
! 69:
! 70: if(MH_Beta != NULL) for (i=0;i<MH_M;i++) b[i]=MH_Beta[i];
! 71: else error_code("MH_Beta is null.");
! 72: if (Ng != NULL) c=a+*(Ng)/2; /* set c, a+n/2, n is given by the option n */
! 73: else error_code("Ng is null.");
! 74:
! 75: bsum = 0; for (i=0; i<m; i++) bsum += b[i];
! 76:
! 77: for (i=0; i<MH_RANK; i++) bitSize[i] = bitcount(i);
! 78: for (i=0, k=0; i<=MH_M; i++) {
! 79: for (j=0; j<MH_RANK; j++) {
! 80: if (bitSize[j] == i) {invbitSize[k] = j; k++; }
! 81: }
! 82: }
! 83: initialized = 1;
! 84: }
! 85: for (i=0; i<MH_M; i++) {
! 86: y[i] = b[i]*x;
! 87: }
! 88: for (i=0; i<MH_M; i++) {
! 89: for (k=0; k<MH_M; k++) {
! 90: if (i != k) {
! 91: yik[idxM(i,k)] = y[k]/(y[i]-y[k]);
! 92: yik2[idxM(i,k)] = y[i]/((y[i]-y[k])*(y[i]-y[k]));
! 93: }
! 94: }
! 95: }
! 96:
! 97: /*
! 98: generation of f2[i][jj] = p(x) pd{i}^2 \pd{J} F
! 99: Do not forget to devide by y[i]
! 100: */
! 101: for (i=0; i<MH_M; i++) for (j=0; j<MH_RANK; j++) f2[idxRank(i,j)] = NaN;
! 102: for (i=0; i<MH_M; i++) {
! 103: f2[idxRank(i,0)] = -((c-y[i])*f[join(i,0)]-a*f[0]);
! 104: for (k=0; k<MH_M; k++) {
! 105: if (i!=k) {
! 106: f2[idxRank(i,0)] += -0.5*yik[idxM(i,k)]*(f[join(i,0)]-f[join(k,0)]);
! 107: }
! 108: }
! 109: f2[idxRank(i,0)] = f2[idxRank(i,0)]/y[i];
! 110: }
! 111: for (p=1; p<MH_RANK; p++) {
! 112: jj = invbitSize[p];
! 113: for (i=0; i<MH_M; i++) {
! 114: if (member(i,jj)) continue;
! 115: ii = join(i,jj);
! 116: rijj = -( (c-y[i])*f[ii]-a*f[jj] );
! 117: for (k=0; k<MH_M; k++) {
! 118: if (!member(k,ii)) rijj += -(0.5)*yik[idxM(i,k)]*(f[ii]-f[join(k,jj)]);
! 119: if (member(k,jj)) {
! 120: rijj += -(0.5)*yik[idxM(i,k)]*f[ii];
! 121: rijj += -(0.5)*yik2[idxM(i,k)]*(f[join(i,delete(k,jj))]-f[jj]);
! 122: }
! 123: }
! 124: f2[idxRank(i,jj)] = rijj;
! 125: for (k = 0; k<MH_M; k++) {
! 126: if (member(k,jj)) {
! 127: f2[idxRank(i,jj)] += (0.5)*yik[idxM(i,k)]*f2[idxRank(k,delete(k,jj))];
! 128: }
! 129: }
! 130: f2[idxRank(i,jj)] = f2[idxRank(i,jj)]/y[i];
! 131: }
! 132: }
! 133: /** showf2("f2",f2,MH_M,MH_RANK); exit(0); */
! 134:
! 135: for (jj=0; jj<MH_RANK; jj++) {
! 136: val[jj] = 0;
! 137: for (i=0; i<MH_M; i++) {
! 138: if (member(i,jj)) {
! 139: val[jj] += b[i]*f2[idxRank(i,delete(i,jj))];
! 140: }else{
! 141: val[jj] += b[i]*f[join(i,jj)];
! 142: }
! 143: }
! 144: }
! 145:
! 146: /* diagonal shift. check MH_M.*/
! 147: for (i=0; i<MH_RANK; i++) {
! 148: val[i] += (-bsum+MH_M*(c-a)/x)*f[i];
! 149: }
! 150: }
! 151:
! 152: error_code(char *s) {
! 153: fprintf(stderr,"%s",s);
! 154: exit(10);
! 155: }
! 156:
! 157: /* for debug */
! 158: showf(char *s,double *v,int m) {
! 159: int i;
! 160: printf("%s=\n",s);
! 161: for (i=0; i<m; i++) {
! 162: printf("%e, ",v[i]);
! 163: }
! 164: printf("\n");
! 165: }
! 166:
! 167: showd(char *s,int *v,int m) {
! 168: int i;
! 169: printf("%s=\n",s);
! 170: for (i=0; i<m; i++) {
! 171: printf("%5d, ",v[i]);
! 172: }
! 173: printf("\n");
! 174: }
! 175:
! 176: showf2(char *s,double *v,int m,int n) {
! 177: int i,j;
! 178: printf("%s=\n",s);
! 179: for (i=0; i<m; i++) {
! 180: for (j=0; j<n; j++) {
! 181: printf("%e, ",v[i*n+j]);
! 182: }
! 183: printf("\n");
! 184: }
! 185: printf("\n");
! 186: printf("member(0,0)=%d, member(0,1)=%d, member(0,2)=%d,member(0,3)=%d\n",
! 187: member(0,0),member(0,1),member(0,2),member(0,3));
! 188: }
! 189:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to code-n.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