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