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

Diff for /OpenXM/src/hgm/mh/src/jack-n.c between version 1.24 and 1.26

version 1.24, 2014/03/16 00:00:46 version 1.26, 2014/03/17 02:49:39
Line 5 
Line 5 
 #include <string.h>  #include <string.h>
 #include "sfile.h"  #include "sfile.h"
 /*  /*
   $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.23 2014/03/15 11:23:58 takayama Exp $    $OpenXM: OpenXM/src/hgm/mh/src/jack-n.c,v 1.25 2014/03/16 03:11:07 takayama Exp $
   Ref: copied from this11/misc-2011/A1/wishart/Prog    Ref: copied from this11/misc-2011/A1/wishart/Prog
   jack-n.c, translated from mh.rr or tk_jack.rr in the asir-contrib. License: LGPL    jack-n.c, translated from mh.rr or tk_jack.rr in the asir-contrib. License: LGPL
   Koev-Edelman for higher order derivatives.    Koev-Edelman for higher order derivatives.
Line 100  static double M_rel_error=0.0; /* relative errors */
Line 100  static double M_rel_error=0.0; /* relative errors */
  If automatic == 1, then the series is reevaluated as long as t_success!=1   If automatic == 1, then the series is reevaluated as long as t_success!=1
  by increasing X0g (evaluation point) and M_m (approx degree);   by increasing X0g (evaluation point) and M_m (approx degree);
  */   */
 static int M_automatic=0;  static int M_automatic=1;
 /* Estimated degree bound for series expansion. See mh_t */  /* Estimated degree bound for series expansion. See mh_t */
 static int M_m_estimated_approx_deg=0;  static int M_m_estimated_approx_deg=0;
 /* Let F(i) be the approximation up to degree i.  /* Let F(i) be the approximation up to degree i.
Line 112  static double M_series_error;
Line 112  static double M_series_error;
   M_series_error < M_assigend_series_error (A) is required for the    M_series_error < M_assigend_series_error (A) is required for the
   estimated_approx_deg.    estimated_approx_deg.
  */   */
 static double M_assigned_series_error=0.00001;  static double M_assigned_series_error=M_ASSIGNED_SERIES_ERROR_DEFAULT;
 /*  /*
   Let Ef be the exponential factor ( Ef=(4)/1F1 of [HNTT] )    Let Ef be the exponential factor ( Ef=(4)/1F1 of [HNTT] )
   If F(M_m)*Ef < x0value_min (B), the success=0 and X0g is increased.    If F(M_m)*Ef < x0value_min (B), the success=0 and X0g is increased.
   Note that minimal double is about 2e-308    Note that minimal double is about 2e-308
  */   */
 static double M_x0value_min=1e-30;  static double M_x0value_min=1e-60;
 /*  /*
   estimated_X0g is the suggested value of X0g.    estimated_X0g is the suggested value of X0g.
  */   */
Line 149  static double M_beta_i_beta_j_min;
Line 149  static double M_beta_i_beta_j_min;
   Value of matrix hg    Value of matrix hg
 */  */
 static double M_mh_t_value;  static double M_mh_t_value;
   /*
    Show the process of updating degree.
    */
   int M_show_autosteps=1;
   
   
 /* prototypes */  /* prototypes */
 static void *mymalloc(int size);  static void *mymalloc(int size);
 static int myfree(void *p);  static int myfree(void *p);
Line 1377  double mh_t(double A[A_LEN],double B[B_LEN],int N,int 
Line 1380  double mh_t(double A[A_LEN],double B[B_LEN],int N,int 
   
   M_series_error = serror;    M_series_error = serror;
   M_recommended_abserr = iv*M_assigned_series_error;    M_recommended_abserr = iv*M_assigned_series_error;
   
   printf("%%%%serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m);  
   printf("%%%%x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);  
   printf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);  
   fprintf(stderr,"%%%%(stderr) serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m);  
   fprintf(stderr,"%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);  
   fprintf(stderr,"%%%%(stderr) F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);  
   
     if (M_show_autosteps) {
       printf("%%%%serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m);
       printf("%%%%x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
       printf("%%%%F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
       fprintf(stderr,"%%%%(stderr) serror=%lg, M_assigned_series_error=%lg, M_m_estimated_approx_deg=%d,M_m=%d\n",serror,M_assigned_series_error,M_m_estimated_approx_deg,M_m);
       fprintf(stderr,"%%%%(stderr) x0value_min=%lg, x0g_bound=%lg\n",M_x0value_min, M_X0g_bound);
       fprintf(stderr,"%%%%(stderr) F=%lg,Ef=%lg,M_estimated_X0g=%lg, X0g=%lg\n",F,iv_factor(),M_estimated_X0g,X0g);
     }
   
   M_mh_t_value=F;    M_mh_t_value=F;
   return(F);    return(F);
 }  }
Line 1544  struct MH_RESULT *jk_main2(int argc,char *argv[],int a
Line 1549  struct MH_RESULT *jk_main2(int argc,char *argv[],int a
     ans->recommended_abserr = 1.0e-10;      ans->recommended_abserr = 1.0e-10;
   }    }
   else ans = NULL;    else ans = NULL;
     if (M_automatic) {
       /* Differentiation can be M_m in the bit pattern in the M_n variable case.*/
       if (M_n > Mapprox) Mapprox=M_n;
     }
   /* Output by a file=stdout */    /* Output by a file=stdout */
   ofp = mh_fopen("stdout","w",JK_byFile);    ofp = mh_fopen("stdout","w",JK_byFile);
   
Line 1630  static int usage() {
Line 1639  static int usage() {
   fprintf(stderr," With the --notable option, it does not use the Lemma 3.2 of Koev-Edelman (there is a typo: kappa'_r = mu'_r for 1<=r<=mu_k).\n");    fprintf(stderr," With the --notable option, it does not use the Lemma 3.2 of Koev-Edelman (there is a typo: kappa'_r = mu'_r for 1<=r<=mu_k).\n");
   fprintf(stderr," An example format of the input_data_file can be obtained by executing hgm_jack-n with no option.\n");    fprintf(stderr," An example format of the input_data_file can be obtained by executing hgm_jack-n with no option.\n");
   fprintf(stderr,"By --automatic option, X0g and degree are automatically searched. The current strategy is described in mh_t in jack-n.c\n");    fprintf(stderr,"By --automatic option, X0g and degree are automatically searched. The current strategy is described in mh_t in jack-n.c\n");
   fprintf(stderr,"Default values for the papameters of the automatic mode: assigned_series_error=%lg, x0value_min=%lg\n",M_assigned_series_error,M_x0value_min);    fprintf(stderr,"Default values for the papameters of the automatic mode: automatic=%d, assigned_series_error=%lg, x0value_min=%lg\n",M_automatic,M_assigned_series_error,M_x0value_min);
   fprintf(stderr,"Todo: automatic mode throws away the table of Jack polynomials of the previous degrees and reevaluate them. They should be kept.\n");    fprintf(stderr,"Todo: automatic mode throws away the table of Jack polynomials of the previous degrees and reevaluate them. They should be kept.\n");
   fprintf(stderr,"\nExamples:\n");    fprintf(stderr,"\nExamples:\n");
   fprintf(stderr,"[1] ./hgm_jack-n \n");    fprintf(stderr,"[1] ./hgm_jack-n \n");
   fprintf(stderr,"[2] ./hgm_jack-n --x0 0.1 \n");    fprintf(stderr,"[2] ./hgm_jack-n --x0 0.1 \n");
   fprintf(stderr,"[3] ./hgm_jack-n --x0 0.1 --degree 15 \n");    fprintf(stderr,"[3] ./hgm_jack-n --x0 0.1 --degree 15\n");
   fprintf(stderr,"[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --degree 15 \n");    fprintf(stderr,"[4] ./hgm_jack-n --idata Testdata/tmp-idata3.txt --degree 15  --automatic 0\n");
   fprintf(stderr,"[5] ./hgm_jack-n --degree 15 >test2.txt\n");    fprintf(stderr,"[5] ./hgm_jack-n --degree 15 >test2.txt\n");
   fprintf(stderr,"    ./hgm_w-n --idata test2.txt --gnuplotf test-g\n");    fprintf(stderr,"    ./hgm_w-n --idata test2.txt --gnuplotf test-g\n");
   fprintf(stderr,"    gnuplot -persist <test-g-gp.txt\n");    fprintf(stderr,"    gnuplot -persist <test-g-gp.txt\n");
Line 1797  static int setParam(char *fname) {
Line 1806  static int setParam(char *fname) {
         fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);          fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
       }        }
       M_X0g_bound = tk.dval;        M_X0g_bound = tk.dval;
         continue;
       }
       if (strcmp(s,"show_autosteps")==0) {
         if (mh_getoken(s,SMAX-1,fp).type != MH_TOKEN_EQ) {
           fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
         }
         if ((tk=mh_getoken(s,SMAX-1,fp)).type != MH_TOKEN_INT) {
           fprintf(stderr,"Syntax error at %s\n",s); mh_exit(-1);
         }
         M_show_autosteps = tk.ival;
       continue;        continue;
     }      }
     fprintf(stderr,"Unknown ID at %s\n",s); mh_exit(-1);      fprintf(stderr,"Unknown ID at %s\n",s); mh_exit(-1);
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.26
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>