Annotation of OpenXM/src/hgm/mh/src/code-n-2f1.c, Revision 1.5
1.1 takayama 1: /*
1.5 ! takayama 2: $OpenXM: OpenXM/src/hgm/mh/src/code-n-2f1.c,v 1.4 2016/02/04 06:56:04 takayama Exp $
1.1 takayama 3: License: LGPL
4: Ref: code-n.c, 2016.01.30, 31.
5: */
6: #include <stdio.h>
7: #include <stdlib.h>
8: #include <math.h>
9: #include "sfile.h"
10: #include "mh.h"
11:
12: static void error_code(char *s);
13: #ifdef STANDALONE
14: static void showf(char *s,double *v,int m);
15: static void showf(char *s,double *v,int m);
16: static void showf2(char *s,double *v,int m,int n);
17: #endif
18:
19: extern int MH_M;
20: extern int MH_RANK;
1.4 takayama 21: extern int MH_Ef_type;
1.5 ! takayama 22: extern double *MH_A_pFq;
! 23: extern double *MH_B_pFq;
1.1 takayama 24:
25: static int bitcount(int n) {int c; c=0; while (n) { ++c; n &= (n-1); } return(c);}
26: #define join(k,jj) ( (1 << k) | jj)
27: #define delete(k,jj) ( (~(1 << k)) & jj)
28: #define member(k,ii) ( (1 << k) & ii)
29: #define NaN 9.999e100
30:
31: #define idxM(i,j) ((i)*MH_M+j)
32: #define idxRank(i,j) ((i)*MH_RANK+j)
33:
34:
1.5 ! takayama 35: void mh_rf_ef_type_2(double x, double *f, int rank_not_used, double *val, int n_not_used)
1.1 takayama 36: { extern double *MH_Beta; /* beta = (1/2)*\simga^(-1) */
37: extern double *MH_Ng; /* freedom */
38: extern int MH_Mg; /* number of variables, MH_Mg=MH_M */
39: static int initialized = 0;
40: extern int MH_deallocate;
41: /* static double b[MH_M]; */
42: static double *b;
43: static double bsum = 0;
44: static int m;
45: /*
46: static int bitSize[MH_RANK];
47: static int invbitSize[MH_RANK];
48: */
49: static int *bitSize;
50: static int *invbitSize;
51: int i,j,k,p;
52: int ii,jj; /* stands for I and J */
53: double rijj;
1.4 takayama 54: double dd;
1.1 takayama 55:
56: static double *pp=NULL; /* p(x_i) */
57: static double *qq=NULL; /* q(x_i,x_k) */
58: static double *qq2=NULL; /* q_2(x_i,x_k) */
59: static double *rr=NULL; /* r(x_i) */
60: static double *qq_pd=NULL; /* dx_k q */
61: static double *qq2_pd=NULL; /* dx_k q_2 */
62: /* double f2[MH_M][MH_RANK];*/
63: static double *f2=NULL;
64: static double *y=NULL;
65: static double aaa=NaN;
66: static double bbb,ccc;
1.4 takayama 67: static double *b2=NULL; /* b_i/(x+b_i)^2 */
1.1 takayama 68:
69: mh_check_intr(100);
70: if (MH_deallocate && initialized) {
71: if (b) mh_free(b);
72: if (bitSize) mh_free(bitSize);
73: if (invbitSize) mh_free(invbitSize);
74: if (f2) mh_free(f2);
75: if (pp) mh_free(pp);
76: if (qq) mh_free(qq);
77: if (qq2) mh_free(qq2);
78: if (rr) mh_free(rr);
79: if (qq_pd) mh_free(qq_pd);
80: if (qq2_pd) mh_free(qq2_pd);
1.4 takayama 81: if (b2) mh_free(b2);
1.1 takayama 82: if (y) mh_free(y);
83: b = f2 = y = NULL;
84: pp = qq = qq2 = rr = qq_pd = qq2_pd = NULL;
85: bitSize = invbitSize = NULL;
86: initialized = 0; return;
87: }
88: if (!initialized) {
89: b = (double *)mh_malloc(sizeof(double)*MH_M);
90: bitSize = (int *)mh_malloc(sizeof(int)*MH_RANK);
91: invbitSize = (int *)mh_malloc(sizeof(int)*MH_RANK);
92: y = (double *)mh_malloc(sizeof(double)*MH_RANK);
93:
94: f2 = (double *)mh_malloc(sizeof(double)*MH_M*MH_RANK);
95: pp = (double *)mh_malloc(sizeof(double)*MH_M);
96: qq = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
97: qq2 = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
98: rr = (double *)mh_malloc(sizeof(double)*MH_M);
99: qq_pd = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
100: qq2_pd = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
1.4 takayama 101: b2 = (double *)mh_malloc(sizeof(double)*MH_M);
1.1 takayama 102:
103: m = MH_Mg;
104: if (m != MH_M) error_code("MH_M != m. MH_M is given by -DMH_M...");
1.5 ! takayama 105: aaa = MH_A_pFq[0];
! 106: bbb = MH_A_pFq[1];
! 107: ccc = MH_B_pFq[0];
1.1 takayama 108:
109: if(MH_Beta != NULL) for (i=0;i<MH_M;i++) b[i]=MH_Beta[i];
110: else error_code("MH_Beta is null.");
111:
112: bsum = 0; for (i=0; i<m; i++) bsum += b[i];
113:
114: for (i=0; i<MH_RANK; i++) bitSize[i] = bitcount(i);
115: for (i=0, k=0; i<=MH_M; i++) {
116: for (j=0; j<MH_RANK; j++) {
117: if (bitSize[j] == i) {invbitSize[k] = j; k++; }
118: }
119: }
120: initialized = 1;
121: }
1.4 takayama 122: if (MH_Ef_type == 2) {
123: for (i=0; i<MH_M; i++) {
124: y[i] = x/(b[i]+x);
125: }
126: }else{
127: for (i=0; i<MH_M; i++) {
128: y[i] = b[i]*x;
129: }
1.1 takayama 130: }
131: for (i=0; i<MH_M; i++) {
132: pp[i] = (ccc-(m-1.0)/2.0-(aaa+bbb+1-(m-1.0)/2.0)*y[i])/(y[i]*(1-y[i]));
133: rr[i] = aaa*bbb/(y[i]*(1-y[i]));
134: for (k=0; k<MH_M; k++) {
135: if (i != k) {
1.2 takayama 136: qq2[idxM(i,k)] = 1.0/(2.0*(y[i]-y[k]));
137: qq[idxM(i,k)] = y[k]*(1-y[k])/(2.0*y[i]*(1-y[i])*(y[i]-y[k]));
138: qq2_pd[idxM(i,k)] = 1.0/(2.0*(y[i]-y[k])*(y[i]-y[k]));
139: qq_pd[idxM(i,k)] = 1/(2.0*y[i]*(1-y[i]));
140: qq_pd[idxM(i,k)] *= ((1-2*y[k])*(y[i]-y[k])+y[k]*(1-y[k]))/((y[i]-y[k])*(y[i]-y[k]));
1.1 takayama 141: }
142: }
143: }
144:
145: /*
146: generation of f2[i][jj] = pd{i}^2 \pd{J} F
147: */
148: for (i=0; i<MH_M; i++) for (j=0; j<MH_RANK; j++) f2[idxRank(i,j)] = NaN;
149: /* The case J = jj = \emptyset */
150: for (i=0; i<MH_M; i++) {
1.3 takayama 151: f2[idxRank(i,0)] = -pp[i]*f[join(i,0)]+rr[i]*f[0];
1.1 takayama 152: for (k=0; k<MH_M; k++) {
153: if (i!=k) {
1.3 takayama 154: f2[idxRank(i,0)] -= qq2[idxM(i,k)]*f[join(i,0)]-qq[idxM(i,k)]*f[join(k,0)];
1.1 takayama 155: }
156: }
157: }
158: for (p=1; p<MH_RANK; p++) {
159: jj = invbitSize[p]; /* evaluate from jj of small bitSize */
160: for (i=0; i<MH_M; i++) {
161: if (member(i,jj)) continue;
162: ii = join(i,jj);
163: rijj = -pp[i]*f[ii]+rr[i]*f[jj];
164: for (k=0; k<MH_M; k++) {
165: if (k==i) continue;
166: rijj -= qq2[idxM(i,k)]*f[ii];
167: if (member(k,jj)) {
1.2 takayama 168: rijj -= qq2_pd[idxM(i,k)]*f[delete(k,ii)]
1.1 takayama 169: -qq_pd[idxM(i,k)]*f[jj];
170: }else{
171: rijj -= -qq[idxM(i,k)]*f[join(k,jj)];
172: }
173: }
174: /* term of the form dx_k^2 dx_{jj-k} */
175: for (k = 0; k<MH_M; k++) {
176: if (k == i) continue;
177: if (member(k,jj)) {
1.2 takayama 178: rijj -= -qq[idxM(i,k)]*f2[idxRank(k,delete(k,jj))];
1.1 takayama 179: }
180: }
181: f2[idxRank(i,jj)] = rijj;
182: }
183: }
184: /** showf2("f2",f2,MH_M,MH_RANK); exit(0); */
185:
1.4 takayama 186: /* sum_j b_j dx_j Base */
187: if (MH_Ef_type==2) {
188: for (i=0; i<MH_M; i++) {
189: b2[i] = b[i]/((b[i]+x)*(b[i]+x));
190: }
191: }else{
192: b2[i] = b[i];
193: }
1.1 takayama 194: for (jj=0; jj<MH_RANK; jj++) {
195: val[jj] = 0;
196: for (i=0; i<MH_M; i++) {
197: if (member(i,jj)) {
1.4 takayama 198: val[jj] += b2[i]*f2[idxRank(i,delete(i,jj))];
1.1 takayama 199: }else{
1.4 takayama 200: val[jj] += b2[i]*f[join(i,jj)];
1.1 takayama 201: }
202: }
203: }
204:
1.4 takayama 205: /* If there is a diagonal shift, add a code .*/
206: if (MH_Ef_type==2) {
207: dd = (ccc-aaa)*(2*aaa-1)/x;
208: for (i=0; i<MH_M; i++) {
209: dd += -bbb/(b[i]+x);
210: }
211: for (jj=0; jj<MH_RANK; jj++) {
212: val[jj] += dd*f[jj];
213: }
214: }
1.1 takayama 215: }
216:
217: static void error_code(char *s) {
218: oxprintfe("%s",s);
219: mh_exit(10);
220: }
221:
222:
223: #ifdef STANDALONE
224: /* for debug */
225: static void showf(char *s,double *v,int m) {
226: int i;
227: oxprintf("%s=\n",s);
228: for (i=0; i<m; i++) {
229: oxprintf("%e, ",v[i]);
230: }
231: oxprintf("\n");
232: }
233:
234: static void showd(char *s,int *v,int m) {
235: int i;
236: oxprintf("%s=\n",s);
237: for (i=0; i<m; i++) {
238: oxprintf("%5d, ",v[i]);
239: }
240: oxprintf("\n");
241: }
242:
243: static void showf2(char *s,double *v,int m,int n) {
244: int i,j;
245: oxprintf("%s=\n",s);
246: for (i=0; i<m; i++) {
247: for (j=0; j<n; j++) {
248: oxprintf("%e, ",v[i*n+j]);
249: }
250: oxprintf("\n");
251: }
252: oxprintf("\n");
253: oxprintf("member(0,0)=%d, member(0,1)=%d, member(0,2)=%d,member(0,3)=%d\n",
254: member(0,0),member(0,1),member(0,2),member(0,3));
255: }
256: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to code-n-2f1.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