Annotation of OpenXM/src/hgm/mh/src/code-n.c, Revision 1.6
1.1 takayama 1: /*
1.6 ! takayama 2: $OpenXM: OpenXM/src/hgm/mh/src/code-n.c,v 1.5 2013/03/05 05:26:07 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.5 takayama 64: mh_check_intr(100);
1.3 takayama 65: if (MH_deallocate && initialized) {
1.6 ! takayama 66: if (b) mh_free(b);
! 67: if (bitSize) mh_free(bitSize);
! 68: if (invbitSize) mh_free(invbitSize);
! 69: if (f2) mh_free(f2);
! 70: if (yik) mh_free(yik);
! 71: if (yik2) mh_free(yik2);
! 72: if (y) mh_free(y);
1.3 takayama 73: b = f2 = yik = yik2 = y = NULL;
74: bitSize = invbitSize = NULL;
1.6 ! takayama 75: initialized = 0; return;
1.3 takayama 76: }
1.1 takayama 77: if (!initialized) {
1.2 takayama 78: b = (double *)mh_malloc(sizeof(double)*MH_M);
79: bitSize = (int *)mh_malloc(sizeof(int)*MH_RANK);
80: invbitSize = (int *)mh_malloc(sizeof(int)*MH_RANK);
81:
82: f2 = (double *)mh_malloc(sizeof(double)*MH_M*MH_RANK);
83: yik = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
84: yik2 = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
85: y = (double *)mh_malloc(sizeof(double)*MH_M);
1.1 takayama 86:
87: m = MH_Mg;
88: if (m != MH_M) error_code("MH_M != m. MH_M is given by -DMH_M...");
89: a = ((double)m+1.0)/2.0;
90:
91: if(MH_Beta != NULL) for (i=0;i<MH_M;i++) b[i]=MH_Beta[i];
92: else error_code("MH_Beta is null.");
1.2 takayama 93: if (MH_Ng != NULL) c=a+*(MH_Ng)/2; /* set c, a+n/2, n is given by the option n */
94: else error_code("MH_Ng is null.");
1.1 takayama 95:
96: bsum = 0; for (i=0; i<m; i++) bsum += b[i];
97:
98: for (i=0; i<MH_RANK; i++) bitSize[i] = bitcount(i);
99: for (i=0, k=0; i<=MH_M; i++) {
100: for (j=0; j<MH_RANK; j++) {
1.6 ! takayama 101: if (bitSize[j] == i) {invbitSize[k] = j; k++; }
1.1 takayama 102: }
103: }
104: initialized = 1;
105: }
106: for (i=0; i<MH_M; i++) {
107: y[i] = b[i]*x;
108: }
109: for (i=0; i<MH_M; i++) {
110: for (k=0; k<MH_M; k++) {
111: if (i != k) {
1.6 ! takayama 112: yik[idxM(i,k)] = y[k]/(y[i]-y[k]);
! 113: yik2[idxM(i,k)] = y[i]/((y[i]-y[k])*(y[i]-y[k]));
1.1 takayama 114: }
115: }
116: }
117:
118: /*
119: generation of f2[i][jj] = p(x) pd{i}^2 \pd{J} F
120: Do not forget to devide by y[i]
121: */
122: for (i=0; i<MH_M; i++) for (j=0; j<MH_RANK; j++) f2[idxRank(i,j)] = NaN;
123: for (i=0; i<MH_M; i++) {
124: f2[idxRank(i,0)] = -((c-y[i])*f[join(i,0)]-a*f[0]);
125: for (k=0; k<MH_M; k++) {
126: if (i!=k) {
1.6 ! takayama 127: f2[idxRank(i,0)] += -0.5*yik[idxM(i,k)]*(f[join(i,0)]-f[join(k,0)]);
1.1 takayama 128: }
129: }
130: f2[idxRank(i,0)] = f2[idxRank(i,0)]/y[i];
131: }
132: for (p=1; p<MH_RANK; p++) {
133: jj = invbitSize[p];
134: for (i=0; i<MH_M; i++) {
135: if (member(i,jj)) continue;
136: ii = join(i,jj);
137: rijj = -( (c-y[i])*f[ii]-a*f[jj] );
138: for (k=0; k<MH_M; k++) {
1.6 ! takayama 139: if (!member(k,ii)) rijj += -(0.5)*yik[idxM(i,k)]*(f[ii]-f[join(k,jj)]);
! 140: if (member(k,jj)) {
! 141: rijj += -(0.5)*yik[idxM(i,k)]*f[ii];
! 142: rijj += -(0.5)*yik2[idxM(i,k)]*(f[join(i,delete(k,jj))]-f[jj]);
! 143: }
1.1 takayama 144: }
145: f2[idxRank(i,jj)] = rijj;
146: for (k = 0; k<MH_M; k++) {
1.6 ! takayama 147: if (member(k,jj)) {
! 148: f2[idxRank(i,jj)] += (0.5)*yik[idxM(i,k)]*f2[idxRank(k,delete(k,jj))];
! 149: }
1.1 takayama 150: }
151: f2[idxRank(i,jj)] = f2[idxRank(i,jj)]/y[i];
152: }
153: }
154: /** showf2("f2",f2,MH_M,MH_RANK); exit(0); */
155:
156: for (jj=0; jj<MH_RANK; jj++) {
157: val[jj] = 0;
158: for (i=0; i<MH_M; i++) {
159: if (member(i,jj)) {
1.6 ! takayama 160: val[jj] += b[i]*f2[idxRank(i,delete(i,jj))];
1.1 takayama 161: }else{
1.6 ! takayama 162: val[jj] += b[i]*f[join(i,jj)];
1.1 takayama 163: }
164: }
165: }
166:
167: /* diagonal shift. check MH_M.*/
168: for (i=0; i<MH_RANK; i++) {
169: val[i] += (-bsum+MH_M*(c-a)/x)*f[i];
170: }
171: }
172:
1.2 takayama 173: static void error_code(char *s) {
1.1 takayama 174: fprintf(stderr,"%s",s);
1.2 takayama 175: mh_exit(10);
1.1 takayama 176: }
177:
178: /* for debug */
1.2 takayama 179: static void showf(char *s,double *v,int m) {
1.1 takayama 180: int i;
181: printf("%s=\n",s);
182: for (i=0; i<m; i++) {
183: printf("%e, ",v[i]);
184: }
185: printf("\n");
186: }
187:
1.2 takayama 188: static void showd(char *s,int *v,int m) {
1.1 takayama 189: int i;
190: printf("%s=\n",s);
191: for (i=0; i<m; i++) {
192: printf("%5d, ",v[i]);
193: }
194: printf("\n");
195: }
196:
197: showf2(char *s,double *v,int m,int n) {
198: int i,j;
199: printf("%s=\n",s);
200: for (i=0; i<m; i++) {
201: for (j=0; j<n; j++) {
202: printf("%e, ",v[i*n+j]);
203: }
204: printf("\n");
205: }
206: printf("\n");
207: printf("member(0,0)=%d, member(0,1)=%d, member(0,2)=%d,member(0,3)=%d\n",
1.6 ! takayama 208: member(0,0),member(0,1),member(0,2),member(0,3));
1.1 takayama 209: }
210:
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