===================================================================
RCS file: /home/cvs/OpenXM/src/hgm/mh/src/code-n-2f1.c,v
retrieving revision 1.2
retrieving revision 1.5
diff -u -p -r1.2 -r1.5
--- OpenXM/src/hgm/mh/src/code-n-2f1.c	2016/02/02 03:00:08	1.2
+++ OpenXM/src/hgm/mh/src/code-n-2f1.c	2016/02/13 02:18:59	1.5
@@ -1,5 +1,5 @@
 /*
-$OpenXM: OpenXM/src/hgm/mh/src/code-n-2f1.c,v 1.1 2016/01/31 02:06:16 takayama Exp $
+$OpenXM: OpenXM/src/hgm/mh/src/code-n-2f1.c,v 1.4 2016/02/04 06:56:04 takayama Exp $
 License: LGPL
 Ref: code-n.c, 2016.01.30, 31.
  */
@@ -9,8 +9,6 @@ Ref: code-n.c, 2016.01.30, 31.
 #include "sfile.h"
 #include "mh.h"
 
-void mh_setabc(double a,double b, double c);
-
 static void error_code(char *s);
 #ifdef STANDALONE
 static void showf(char *s,double *v,int m);
@@ -20,6 +18,9 @@ static void showf2(char *s,double *v,int m,int n);
 
 extern int MH_M;
 extern int MH_RANK;
+extern int MH_Ef_type;
+extern double *MH_A_pFq;
+extern double *MH_B_pFq;
 
 static int bitcount(int n) {int c; c=0; while (n) { ++c; n &= (n-1); } return(c);}
 #define join(k,jj) ( (1 << k) | jj)
@@ -30,10 +31,8 @@ static int bitcount(int n) {int c; c=0; while (n) { ++
 #define idxM(i,j) ((i)*MH_M+j)
 #define idxRank(i,j) ((i)*MH_RANK+j)
 
-static double MH_aaa=NaN;
-static double MH_bbb, MH_ccc;
 
-void mh_rf(double x, double *f, int rank_not_used, double *val, int n_not_used)
+void mh_rf_ef_type_2(double x, double *f, int rank_not_used, double *val, int n_not_used)
 { extern double *MH_Beta; /* beta = (1/2)*\simga^(-1) */
   extern double *MH_Ng;   /* freedom */
   extern int MH_Mg;       /* number of variables, MH_Mg=MH_M */
@@ -52,6 +51,7 @@ void mh_rf(double x, double *f, int rank_not_used, dou
   int i,j,k,p;
   int ii,jj; /* stands for I and J */
   double rijj;
+  double dd;
 
   static double *pp=NULL;  /* p(x_i) */
   static double *qq=NULL;  /* q(x_i,x_k) */
@@ -64,6 +64,7 @@ void mh_rf(double x, double *f, int rank_not_used, dou
   static double *y=NULL;
   static double aaa=NaN;
   static double bbb,ccc;
+  static double *b2=NULL;  /* b_i/(x+b_i)^2 */
 
   mh_check_intr(100);
   if (MH_deallocate && initialized) {
@@ -77,12 +78,12 @@ void mh_rf(double x, double *f, int rank_not_used, dou
         if (rr) mh_free(rr);
         if (qq_pd) mh_free(qq_pd);
         if (qq2_pd) mh_free(qq2_pd);
+		if (b2) mh_free(b2);
 	if (y) mh_free(y);
     b = f2 =  y = NULL;
     pp = qq = qq2 = rr = qq_pd = qq2_pd = NULL;
     bitSize = invbitSize = NULL;
         initialized = 0; return;
-	MH_aaa = NaN;
   }
   if (!initialized) {
     b = (double *)mh_malloc(sizeof(double)*MH_M);
@@ -97,14 +98,13 @@ void mh_rf(double x, double *f, int rank_not_used, dou
     rr   = (double *)mh_malloc(sizeof(double)*MH_M);
     qq_pd = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
     qq2_pd = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
+    b2 = (double *)mh_malloc(sizeof(double)*MH_M);
 
-
     m = MH_Mg;
     if (m != MH_M) error_code("MH_M != m. MH_M is given by -DMH_M...");
-    if (MH_aaa == NaN) error_code("aaa, bbb, ccc are not initialized.");
-    aaa = MH_aaa;
-    bbb = MH_bbb;
-    ccc = MH_ccc;
+    aaa = MH_A_pFq[0];
+    bbb = MH_A_pFq[1];
+    ccc = MH_B_pFq[0];
 
     if(MH_Beta != NULL) for (i=0;i<MH_M;i++) b[i]=MH_Beta[i];
     else error_code("MH_Beta is null.");
@@ -119,8 +119,14 @@ void mh_rf(double x, double *f, int rank_not_used, dou
     }
     initialized = 1;
   }
-  for (i=0; i<MH_M; i++) {
-    y[i] = b[i]*x;    
+  if (MH_Ef_type == 2) {
+    for (i=0; i<MH_M; i++) {
+      y[i] = x/(b[i]+x);    
+    }
+  }else{
+    for (i=0; i<MH_M; i++) {
+      y[i] = b[i]*x;    
+    }
   }
   for (i=0; i<MH_M; i++) {
     pp[i] = (ccc-(m-1.0)/2.0-(aaa+bbb+1-(m-1.0)/2.0)*y[i])/(y[i]*(1-y[i]));
@@ -142,10 +148,10 @@ void mh_rf(double x, double *f, int rank_not_used, dou
   for (i=0; i<MH_M; i++) for (j=0; j<MH_RANK; j++) f2[idxRank(i,j)] = NaN;
   /* The case J = jj = \emptyset */  
   for (i=0; i<MH_M; i++) {
-    f2[idxRank(i,0)] = -pp[i]*f[join(i,0)]-qq2[idxM(i,k)]*f[join(i,0)]+rr[i]*f[0];
+    f2[idxRank(i,0)] = -pp[i]*f[join(i,0)]+rr[i]*f[0];
     for (k=0; k<MH_M; k++) {
       if (i!=k) {
-        f2[idxRank(i,0)] -= -qq[idxM(i,k)]*f[join(k,0)];
+        f2[idxRank(i,0)] -= qq2[idxM(i,k)]*f[join(i,0)]-qq[idxM(i,k)]*f[join(k,0)];
       }
     }
   }
@@ -177,19 +183,35 @@ void mh_rf(double x, double *f, int rank_not_used, dou
   }
   /** showf2("f2",f2,MH_M,MH_RANK); exit(0);  */
 
-  /* sum_j b_j dx_j Base */ 
+  /* sum_j b_j dx_j Base */
+  if (MH_Ef_type==2) {
+	for (i=0; i<MH_M; i++) {
+	  b2[i] = b[i]/((b[i]+x)*(b[i]+x));
+	}
+  }else{
+	b2[i] = b[i];
+  }
   for (jj=0; jj<MH_RANK; jj++) {
     val[jj] = 0;
     for (i=0; i<MH_M; i++) {
       if (member(i,jj)) {
-        val[jj] += b[i]*f2[idxRank(i,delete(i,jj))];
+        val[jj] += b2[i]*f2[idxRank(i,delete(i,jj))];
       }else{
-        val[jj] += b[i]*f[join(i,jj)];
+        val[jj] += b2[i]*f[join(i,jj)];
       }
     }
   }
 
-   /*   If there is a diagonal shift, add a code .*/
+  /*   If there is a diagonal shift, add a code .*/
+  if (MH_Ef_type==2) {
+	dd = (ccc-aaa)*(2*aaa-1)/x;
+    for (i=0; i<MH_M; i++) {
+      dd += -bbb/(b[i]+x);
+	}
+	for (jj=0; jj<MH_RANK; jj++) {
+	  val[jj] += dd*f[jj];
+	}
+  }
 }
 
 static void error_code(char *s) {
@@ -197,12 +219,6 @@ static void error_code(char *s) {
   mh_exit(10);
 }
 
-void mh_setabc(double a,double b, double c) {
-  MH_aaa = a;
-  MH_bbb = b;
-  MH_ccc = c;
-  return;
-}
 
 #ifdef STANDALONE
 /* for debug */