[BACK]Return to code-n-2f1.c CVS log [TXT][DIR] Up to [local] / OpenXM / src / hgm / mh / src

Diff for /OpenXM/src/hgm/mh/src/code-n-2f1.c between version 1.3 and 1.5

version 1.3, 2016/02/02 10:58:23 version 1.5, 2016/02/13 02:18:59
Line 1 
Line 1 
 /*  /*
 $OpenXM: OpenXM/src/hgm/mh/src/code-n-2f1.c,v 1.2 2016/02/02 03:00:08 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  License: LGPL
 Ref: code-n.c, 2016.01.30, 31.  Ref: code-n.c, 2016.01.30, 31.
  */   */
Line 9  Ref: code-n.c, 2016.01.30, 31.
Line 9  Ref: code-n.c, 2016.01.30, 31.
 #include "sfile.h"  #include "sfile.h"
 #include "mh.h"  #include "mh.h"
   
 void mh_setabc(double a,double b, double c);  
   
 static void error_code(char *s);  static void error_code(char *s);
 #ifdef STANDALONE  #ifdef STANDALONE
 static void showf(char *s,double *v,int m);  static void showf(char *s,double *v,int m);
Line 20  static void showf2(char *s,double *v,int m,int n);
Line 18  static void showf2(char *s,double *v,int m,int n);
   
 extern int MH_M;  extern int MH_M;
 extern int MH_RANK;  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);}  static int bitcount(int n) {int c; c=0; while (n) { ++c; n &= (n-1); } return(c);}
 #define join(k,jj) ( (1 << k) | jj)  #define join(k,jj) ( (1 << k) | jj)
Line 30  static int bitcount(int n) {int c; c=0; while (n) { ++
Line 31  static int bitcount(int n) {int c; c=0; while (n) { ++
 #define idxM(i,j) ((i)*MH_M+j)  #define idxM(i,j) ((i)*MH_M+j)
 #define idxRank(i,j) ((i)*MH_RANK+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_Beta; /* beta = (1/2)*\simga^(-1) */
   extern double *MH_Ng;   /* freedom */    extern double *MH_Ng;   /* freedom */
   extern int MH_Mg;       /* number of variables, MH_Mg=MH_M */    extern int MH_Mg;       /* number of variables, MH_Mg=MH_M */
Line 52  void mh_rf(double x, double *f, int rank_not_used, dou
Line 51  void mh_rf(double x, double *f, int rank_not_used, dou
   int i,j,k,p;    int i,j,k,p;
   int ii,jj; /* stands for I and J */    int ii,jj; /* stands for I and J */
   double rijj;    double rijj;
     double dd;
   
   static double *pp=NULL;  /* p(x_i) */    static double *pp=NULL;  /* p(x_i) */
   static double *qq=NULL;  /* q(x_i,x_k) */    static double *qq=NULL;  /* q(x_i,x_k) */
Line 64  void mh_rf(double x, double *f, int rank_not_used, dou
Line 64  void mh_rf(double x, double *f, int rank_not_used, dou
   static double *y=NULL;    static double *y=NULL;
   static double aaa=NaN;    static double aaa=NaN;
   static double bbb,ccc;    static double bbb,ccc;
     static double *b2=NULL;  /* b_i/(x+b_i)^2 */
   
   mh_check_intr(100);    mh_check_intr(100);
   if (MH_deallocate && initialized) {    if (MH_deallocate && initialized) {
Line 77  void mh_rf(double x, double *f, int rank_not_used, dou
Line 78  void mh_rf(double x, double *f, int rank_not_used, dou
         if (rr) mh_free(rr);          if (rr) mh_free(rr);
         if (qq_pd) mh_free(qq_pd);          if (qq_pd) mh_free(qq_pd);
         if (qq2_pd) mh_free(qq2_pd);          if (qq2_pd) mh_free(qq2_pd);
                   if (b2) mh_free(b2);
         if (y) mh_free(y);          if (y) mh_free(y);
     b = f2 =  y = NULL;      b = f2 =  y = NULL;
     pp = qq = qq2 = rr = qq_pd = qq2_pd = NULL;      pp = qq = qq2 = rr = qq_pd = qq2_pd = NULL;
     bitSize = invbitSize = NULL;      bitSize = invbitSize = NULL;
         initialized = 0; return;          initialized = 0; return;
         MH_aaa = NaN;  
   }    }
   if (!initialized) {    if (!initialized) {
     b = (double *)mh_malloc(sizeof(double)*MH_M);      b = (double *)mh_malloc(sizeof(double)*MH_M);
Line 97  void mh_rf(double x, double *f, int rank_not_used, dou
Line 98  void mh_rf(double x, double *f, int rank_not_used, dou
     rr   = (double *)mh_malloc(sizeof(double)*MH_M);      rr   = (double *)mh_malloc(sizeof(double)*MH_M);
     qq_pd = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);      qq_pd = (double *)mh_malloc(sizeof(double)*MH_M*MH_M);
     qq2_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;      m = MH_Mg;
     if (m != MH_M) error_code("MH_M != m. MH_M is given by -DMH_M...");      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_A_pFq[0];
     aaa = MH_aaa;      bbb = MH_A_pFq[1];
     bbb = MH_bbb;      ccc = MH_B_pFq[0];
     ccc = MH_ccc;  
   
     if(MH_Beta != NULL) for (i=0;i<MH_M;i++) b[i]=MH_Beta[i];      if(MH_Beta != NULL) for (i=0;i<MH_M;i++) b[i]=MH_Beta[i];
     else error_code("MH_Beta is null.");      else error_code("MH_Beta is null.");
Line 119  void mh_rf(double x, double *f, int rank_not_used, dou
Line 119  void mh_rf(double x, double *f, int rank_not_used, dou
     }      }
     initialized = 1;      initialized = 1;
   }    }
   for (i=0; i<MH_M; i++) {    if (MH_Ef_type == 2) {
     y[i] = b[i]*x;      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++) {    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]));      pp[i] = (ccc-(m-1.0)/2.0-(aaa+bbb+1-(m-1.0)/2.0)*y[i])/(y[i]*(1-y[i]));
Line 177  void mh_rf(double x, double *f, int rank_not_used, dou
Line 183  void mh_rf(double x, double *f, int rank_not_used, dou
   }    }
   /** showf2("f2",f2,MH_M,MH_RANK); exit(0);  */    /** 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++) {    for (jj=0; jj<MH_RANK; jj++) {
     val[jj] = 0;      val[jj] = 0;
     for (i=0; i<MH_M; i++) {      for (i=0; i<MH_M; i++) {
       if (member(i,jj)) {        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{        }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) {  static void error_code(char *s) {
Line 197  static void error_code(char *s) {
Line 219  static void error_code(char *s) {
   mh_exit(10);    mh_exit(10);
 }  }
   
 void mh_setabc(double a,double b, double c) {  
   MH_aaa = a;  
   MH_bbb = b;  
   MH_ccc = c;  
   return;  
 }  
   
 #ifdef STANDALONE  #ifdef STANDALONE
 /* for debug */  /* for debug */
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.5
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>