[BACK]Return to imach.c CVS log [TXT] [DIR] Up to [Local Repository] / imach / src

Diff for /imach/src/imach.c between versions 1.175 and 1.184

version 1.175, 2015/01/03 16:33:42 version 1.184, 2015/03/11 11:52:39
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.184  2015/03/11 11:52:39  brouard
     Summary: Back from Windows 8. Intel Compiler
   
     Revision 1.183  2015/03/10 20:34:32  brouard
     Summary: 0.98q0, trying with directest, mnbrak fixed
   
     We use directest instead of original Powell test; probably no
     incidence on the results, but better justifications;
     We fixed Numerical Recipes mnbrak routine which was wrong and gave
     wrong results.
   
     Revision 1.182  2015/02/12 08:19:57  brouard
     Summary: Trying to keep directest which seems simpler and more general
     Author: Nicolas Brouard
   
     Revision 1.181  2015/02/11 23:22:24  brouard
     Summary: Comments on Powell added
   
     Author:
   
     Revision 1.180  2015/02/11 17:33:45  brouard
     Summary: Finishing move from main to function (hpijx and prevalence_limit)
   
     Revision 1.179  2015/01/04 09:57:06  brouard
     Summary: back to OS/X
   
     Revision 1.178  2015/01/04 09:35:48  brouard
     *** empty log message ***
   
     Revision 1.177  2015/01/03 18:40:56  brouard
     Summary: Still testing ilc32 on OSX
   
     Revision 1.176  2015/01/03 16:45:04  brouard
     *** empty log message ***
   
   Revision 1.175  2015/01/03 16:33:42  brouard    Revision 1.175  2015/01/03 16:33:42  brouard
   *** empty log message ***    *** empty log message ***
   
Line 541 Line 576
 */  */
   
 #define POWELL /* Instead of NLOPT */  #define POWELL /* Instead of NLOPT */
   /* #define POWELLORIGINAL */ /* Don't use Directest to decide new direction but original Powell test */
   /* #define MNBRAKORIGINAL */ /* Don't use mnbrak fix */
   
 #include <math.h>  #include <math.h>
 #include <stdio.h>  #include <stdio.h>
Line 626  typedef struct { Line 663  typedef struct {
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.98p, December 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";  char version[]="Imach version 0.98q0, March 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 742  static double maxarg1,maxarg2; Line 779  static double maxarg1,maxarg2;
 #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))  #define SIGN(a,b) ((b)>0.0 ? fabs(a) : -fabs(a))
 #define rint(a) floor(a+0.5)  #define rint(a) floor(a+0.5)
 /* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */  /* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/myutils_8h-source.html */
 /* #define mytinydouble 1.0e-16 */  #define mytinydouble 1.0e-16
 /* #define DEQUAL(a,b) (fabs((a)-(b))<mytinydouble) */  /* #define DEQUAL(a,b) (fabs((a)-(b))<mytinydouble) */
 /* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */  /* http://www.thphys.uni-heidelberg.de/~robbers/cmbeasy/doc/html/mynrutils_8h-source.html */
 /* static double dsqrarg; */  /* static double dsqrarg; */
Line 801  static int split( char *path, char *dirc Line 838  static int split( char *path, char *dirc
       printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/        printf("Warning you should use %s as a separator\n",DIRSEPARATOR);*/
     /* get current working directory */      /* get current working directory */
     /*    extern  char* getcwd ( char *buf , int len);*/      /*    extern  char* getcwd ( char *buf , int len);*/
     if ( getcwd( dirc, FILENAME_MAX ) == NULL ) {  #ifdef WIN32
       if (_getcwd( dirc, FILENAME_MAX ) == NULL ) {
   #else
           if (getcwd(dirc, FILENAME_MAX) == NULL) {
   #endif
       return( GLOCK_ERROR_GETCWD );        return( GLOCK_ERROR_GETCWD );
     }      }
     /* got dirc from getcwd*/      /* got dirc from getcwd*/
Line 1258  double brent(double ax, double bx, doubl Line 1299  double brent(double ax, double bx, doubl
     if (fu <= fx) {       if (fu <= fx) { 
       if (u >= x) a=x; else b=x;         if (u >= x) a=x; else b=x; 
       SHFT(v,w,x,u)         SHFT(v,w,x,u) 
         SHFT(fv,fw,fx,fu)         SHFT(fv,fw,fx,fu) 
         } else {       } else { 
           if (u < x) a=u; else b=u;         if (u < x) a=u; else b=u; 
           if (fu <= fw || w == x) {         if (fu <= fw || w == x) { 
             v=w;           v=w; 
             w=u;           w=u; 
             fv=fw;           fv=fw; 
             fw=fu;           fw=fu; 
           } else if (fu <= fv || v == x || v == w) {         } else if (fu <= fv || v == x || v == w) { 
             v=u;           v=u; 
             fv=fu;           fv=fu; 
           }         } 
         }       } 
   }     } 
   nrerror("Too many iterations in brent");     nrerror("Too many iterations in brent"); 
   *xmin=x;     *xmin=x; 
Line 1281  double brent(double ax, double bx, doubl Line 1322  double brent(double ax, double bx, doubl
   
 void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,   void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc, 
             double (*func)(double))               double (*func)(double)) 
 {   { /* Given a function func , and given distinct initial points ax and bx , this routine searches in
   the downhill direction (defined by the function as evaluated at the initial points) and returns
   new points ax , bx , cx that bracket a minimum of the function. Also returned are the function
   values at the three points, fa, fb , and fc such that fa > fb and fb < fc.
      */
   double ulim,u,r,q, dum;    double ulim,u,r,q, dum;
   double fu;     double fu; 
     
Line 1289  void mnbrak(double *ax, double *bx, doub Line 1334  void mnbrak(double *ax, double *bx, doub
   *fb=(*func)(*bx);     *fb=(*func)(*bx); 
   if (*fb > *fa) {     if (*fb > *fa) { 
     SHFT(dum,*ax,*bx,dum)       SHFT(dum,*ax,*bx,dum) 
       SHFT(dum,*fb,*fa,dum)       SHFT(dum,*fb,*fa,dum) 
       }     } 
   *cx=(*bx)+GOLD*(*bx-*ax);     *cx=(*bx)+GOLD*(*bx-*ax); 
   *fc=(*func)(*cx);     *fc=(*func)(*cx); 
   while (*fb > *fc) { /* Declining fa, fb, fc */  #ifdef DEBUG
     printf("mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);
     fprintf(ficlog,"mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);
   #endif
     while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc */
     r=(*bx-*ax)*(*fb-*fc);       r=(*bx-*ax)*(*fb-*fc); 
     q=(*bx-*cx)*(*fb-*fa);       q=(*bx-*cx)*(*fb-*fa); 
     u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/       u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ 
       (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscisse of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */        (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); /* Minimum abscissa of a parabolic estimated from (a,fa), (b,fb) and (c,fc). */
     ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscisse where function can be evaluated */      ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscissa where function should be evaluated */
     if ((*bx-u)*(u-*cx) > 0.0) { /* if u between b and c */      if ((*bx-u)*(u-*cx) > 0.0) { /* if u_p is between b and c */
       fu=(*func)(u);         fu=(*func)(u); 
 #ifdef DEBUG  #ifdef DEBUG
       /* f(x)=A(x-u)**2+f(u) */        /* f(x)=A(x-u)**2+f(u) */
Line 1308  void mnbrak(double *ax, double *bx, doub Line 1357  void mnbrak(double *ax, double *bx, doub
       fparabu= *fa - A*(*ax-u)*(*ax-u);        fparabu= *fa - A*(*ax-u)*(*ax-u);
       printf("mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);        printf("mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
       fprintf(ficlog, "mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);        fprintf(ficlog, "mnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu);
         /* And thus,it can be that fu > *fc even if fparabu < *fc */
         /* mnbrak (*ax=7.666299858533, *fa=299039.693133272231), (*bx=8.595447774979, *fb=298976.598289369489),
           (*cx=10.098840694817, *fc=298946.631474258087),  (*u=9.852501168332, fu=298948.773013752128, fparabu=298945.434711494134) */
         /* In that case, there is no bracket in the output! Routine is wrong with many consequences.*/
 #endif   #endif 
   #ifdef MNBRAKORIGINAL
   #else
         if (fu > *fc) {
   #ifdef DEBUG
         printf("mnbrak4  fu > fc \n");
         fprintf(ficlog, "mnbrak4 fu > fc\n");
   #endif
           /* SHFT(u,*cx,*cx,u) /\* ie a=c, c=u and u=c; in that case, next SHFT(a,b,c,u) will give a=b=b, b=c=u, c=u=c and *\/  */
           /* SHFT(*fa,*fc,fu,*fc) /\* (b, u, c) is a bracket while test fb > fc will be fu > fc  will exit *\/ */
           dum=u; /* Shifting c and u */
           u = *cx;
           *cx = dum;
           dum = fu;
           fu = *fc;
           *fc =dum;
         } else { /* end */
   #ifdef DEBUG
         printf("mnbrak3  fu < fc \n");
         fprintf(ficlog, "mnbrak3 fu < fc\n");
   #endif
           dum=u; /* Shifting c and u */
           u = *cx;
           *cx = dum;
           dum = fu;
           fu = *fc;
           *fc =dum;
         }
   #endif
     } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */      } else if ((*cx-u)*(u-ulim) > 0.0) { /* u is after c but before ulim */
   #ifdef DEBUG
         printf("mnbrak2  u after c but before ulim\n");
         fprintf(ficlog, "mnbrak2 u after c but before ulim\n");
   #endif
       fu=(*func)(u);         fu=(*func)(u); 
       if (fu < *fc) {         if (fu < *fc) { 
   #ifdef DEBUG
         printf("mnbrak2  u after c but before ulim AND fu < fc\n");
         fprintf(ficlog, "mnbrak2 u after c but before ulim AND fu <fc \n");
   #endif
         SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))           SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) 
           SHFT(*fb,*fc,fu,(*func)(u))           SHFT(*fb,*fc,fu,(*func)(u)) 
           }         } 
     } else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */      } else if ((u-ulim)*(ulim-*cx) >= 0.0) { /* u outside ulim (verifying that ulim is beyond c) */
   #ifdef DEBUG
         printf("mnbrak2  u outside ulim (verifying that ulim is beyond c)\n");
         fprintf(ficlog, "mnbrak2 u outside ulim (verifying that ulim is beyond c)\n");
   #endif
       u=ulim;         u=ulim; 
       fu=(*func)(u);         fu=(*func)(u); 
     } else {       } else { /* u could be left to b (if r > q parabola has a maximum) */
   #ifdef DEBUG
         printf("mnbrak2  u could be left to b (if r > q parabola has a maximum)\n");
         fprintf(ficlog, "mnbrak2  u could be left to b (if r > q parabola has a maximum)\n");
   #endif
       u=(*cx)+GOLD*(*cx-*bx);         u=(*cx)+GOLD*(*cx-*bx); 
       fu=(*func)(u);         fu=(*func)(u); 
     }       } /* end tests */
     SHFT(*ax,*bx,*cx,u)       SHFT(*ax,*bx,*cx,u) 
       SHFT(*fa,*fb,*fc,fu)       SHFT(*fa,*fb,*fc,fu) 
       }   #ifdef DEBUG
         printf("mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu);
         fprintf(ficlog, "mnbrak2 (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu);
   #endif
     } /* end while; ie return (a, b, c, fa, fb, fc) such that a < b < c with f(a) > f(b) and fb < f(c) */
 }   } 
   
 /*************** linmin ************************/  /*************** linmin ************************/
Line 1389  void powell(double p[], double **xi, int Line 1490  void powell(double p[], double **xi, int
               double (*func)(double []));                 double (*func)(double [])); 
   int i,ibig,j;     int i,ibig,j; 
   double del,t,*pt,*ptt,*xit;    double del,t,*pt,*ptt,*xit;
     double directest;
   double fp,fptt;    double fp,fptt;
   double *xits;    double *xits;
   int niterf, itmp;    int niterf, itmp;
Line 1448  void powell(double p[], double **xi, int Line 1550  void powell(double p[], double **xi, int
 #endif  #endif
       printf("%d",i);fflush(stdout);        printf("%d",i);fflush(stdout);
       fprintf(ficlog,"%d",i);fflush(ficlog);        fprintf(ficlog,"%d",i);fflush(ficlog);
       linmin(p,xit,n,fret,func);         linmin(p,xit,n,fret,func); /* xit[n] has been loaded for direction i */
       if (fabs(fptt-(*fret)) > del) {         if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions 
                                          because that direction will be replaced unless the gain del is small
                                         in comparison with the 'probable' gain, mu^2, with the last average direction.
                                         Unless the n directions are conjugate some gain in the determinant may be obtained
                                         with the new direction.
                                         */
         del=fabs(fptt-(*fret));           del=fabs(fptt-(*fret)); 
         ibig=i;           ibig=i; 
       }         } 
Line 1469  void powell(double p[], double **xi, int Line 1576  void powell(double p[], double **xi, int
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
 #endif  #endif
     } /* end i */      } /* end i */
     if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {      if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */
 #ifdef DEBUG  #ifdef DEBUG
       int k[2],l;        int k[2],l;
       k[0]=1;        k[0]=1;
Line 1501  void powell(double p[], double **xi, int Line 1608  void powell(double p[], double **xi, int
       return;         return; 
     }       } 
     if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");       if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); 
     for (j=1;j<=n;j++) { /* Computes an extrapolated point */      for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
       ptt[j]=2.0*p[j]-pt[j];         ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j];         xit[j]=p[j]-pt[j]; 
       pt[j]=p[j];         pt[j]=p[j]; 
     }       } 
     fptt=(*func)(ptt);       fptt=(*func)(ptt); /* f_3 */
     if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */      if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
       /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */        /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */
       /* From x1 (P0) distance of x2 is at h and x3 is 2h */        /* From x1 (P0) distance of x2 is at h and x3 is 2h */
       /* Let f"(x2) be the 2nd derivative equal everywhere.  */        /* Let f"(x2) be the 2nd derivative equal everywhere.  */
       /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */        /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */
       /* will reach at f3 = fm + h^2/2 f"m  ; f" = (f1 -2f2 +f3 ) / h**2 */        /* will reach at f3 = fm + h^2/2 f"m  ; f" = (f1 -2f2 +f3 ) / h**2 */
       /* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */        /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */
       /* Thus we compare delta(2h) with observed f1-f3 */  
       /* or best gain on one ancient line 'del' with total  */  
       /* gain f1-f2 = f1 - f2 - 'del' with del  */  
       /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */        /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
   #ifdef NRCORIGINAL
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del);        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
   #else
         t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del); /* Intel compiler doesn't work on one line; bug reported */
       t= t- del*SQR(fp-fptt);        t= t- del*SQR(fp-fptt);
       printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);  #endif
       fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);        directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */
 #ifdef DEBUG  #ifdef DEBUG
         printf("t1= %.12lf, t2= %.12lf, t=%.12lf  directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
         fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
       printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),        printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
              (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));               (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
       fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),        fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
Line 1531  void powell(double p[], double **xi, int Line 1639  void powell(double p[], double **xi, int
       printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);        printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
       fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);        fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
 #endif  #endif
       if (t < 0.0) { /* Then we use it for last direction */  #ifdef POWELLORIGINAL
         linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/        if (t < 0.0) { /* Then we use it for new direction */
   #else
         if (directest*t < 0.0) { /* Contradiction between both tests */
         printf("directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);
         printf("f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
         fprintf(ficlog,"directest= %.12lf, t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);
         fprintf(ficlog,"f1-2f2+f3= %.12lf, f1-f2-del= %.12lf, f1-f3= %.12lf\n",fp-2.0*(*fret)+fptt, fp -(*fret) -del, fp-fptt);
       } 
         if (directest < 0.0) { /* Then we use it for new direction */
   #endif
           linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
           xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */            xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
           xi[j][n]=xit[j];      /* and nth direction by the extrapolated */            xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
         }          }
         printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
         fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
   
 #ifdef DEBUG  #ifdef DEBUG
         printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);          printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
Line 1901  double func( double *x) Line 2019  double func( double *x)
         which slows down the processing. The difference can be up to 10%          which slows down the processing. The difference can be up to 10%
         lower mortality.          lower mortality.
           */            */
           lli=log(out[s1][s2] - savm[s1][s2]);          /* If, at the beginning of the maximization mostly, the
              cumulative probability or probability to be dead is
              constant (ie = 1) over time d, the difference is equal to
              0.  out[s1][3] = savm[s1][3]: probability, being at state
              s1 at precedent wave, to be dead a month before current
              wave is equal to probability, being at state s1 at
              precedent wave, to be dead at mont of the current
              wave. Then the observed probability (that this person died)
              is null according to current estimated parameter. In fact,
              it should be very low but not zero otherwise the log go to
              infinity.
           */
   /* #ifdef INFINITYORIGINAL */
   /*          lli=log(out[s1][s2] - savm[s1][s2]); */
   /* #else */
   /*        if ((out[s1][s2] - savm[s1][s2]) < mytinydouble)  */
   /*          lli=log(mytinydouble); */
   /*        else */
   /*          lli=log(out[s1][s2] - savm[s1][s2]); */
   /* #endif */
               lli=log(out[s1][s2] - savm[s1][s2]);
   
         } else if  (s2==-2) {          } else if  (s2==-2) {
           for (j=1,survp=0. ; j<=nlstate; j++)             for (j=1,survp=0. ; j<=nlstate; j++) 
Line 1933  double func( double *x) Line 2070  double func( double *x)
         ipmx +=1;          ipmx +=1;
         sw += weight[i];          sw += weight[i];
         ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;          ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
           /* if (lli < log(mytinydouble)){ */
           /*   printf("Close to inf lli = %.10lf <  %.10lf i= %d mi= %d, s[%d][i]=%d s1=%d s2=%d\n", lli,log(mytinydouble), i, mi,mw[mi][i], s[mw[mi][i]][i], s1,s2); */
           /*   fprintf(ficlog,"Close to inf lli = %.10lf i= %d mi= %d, s[mw[mi][i]][i]=%d\n", lli, i, mi,s[mw[mi][i]][i]); */
           /* } */
       } /* end of wave */        } /* end of wave */
     } /* end of individual */      } /* end of individual */
   }  else if(mle==2){    }  else if(mle==2){
Line 2268  void mlikeli(FILE *ficres,double p[], in Line 2409  void mlikeli(FILE *ficres,double p[], in
 #endif  #endif
   free_matrix(xi,1,npar,1,npar);    free_matrix(xi,1,npar,1,npar);
   fclose(ficrespow);    fclose(ficrespow);
   printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));    printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
   fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));    fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
   fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));    fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
   
 }  }
   
Line 5501  BOOL IsWow64() Line 5642  BOOL IsWow64()
         return bIsWow64 != FALSE;          return bIsWow64 != FALSE;
 }  }
 #endif  #endif
   
 void syscompilerinfo()  void syscompilerinfo()
  {   {
    /* #include "syscompilerinfo.h"*/     /* #include "syscompilerinfo.h"*/
 #if defined(__GNUC__)           /* command line Intel compiler 64bit windows:
           /GS /W3 /Gy /Zc:wchar_t /Zi /O2 /Fd"x64\Release\vc120.pdb" /D "WIN32" 
           /D "NDEBUG" /D "_CONSOLE" /D "_LIB" /D "_UNICODE" /D "UNICODE" /Qipo 
           /Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo 
           /Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch" */
           /*
           /GS /W3 /Gy /Zc:wchar_t /Zi /O3 /Fd"x64\Release\vc120.pdb" /D "WIN32" 
           /D "NDEBUG" /D "_CONSOLE" /D "_LIB" /D "_UNICODE" /D "UNICODE" /Qipo 
           /Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo /Qparallel 
           /Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch" */
   #if defined __INTEL_COMPILER
   #if defined(__GNUC__)
           struct utsname sysInfo;  /* For Intel on Linux and OS/X */
   #endif
   #elif defined(__GNUC__) 
   #ifndef  __APPLE__
 #include <gnu/libc-version.h>  /* Only on gnu */  #include <gnu/libc-version.h>  /* Only on gnu */
 #endif  #endif
      struct utsname sysInfo;
 #include <stdint.h>  
    int cross = CROSS;     int cross = CROSS;
    if (cross){     if (cross){
      printf("Cross-");             printf("Cross-");
      fprintf(ficlog,"Cross-");             fprintf(ficlog, "Cross-");
    }     }
   #endif
   
   #include <stdint.h>
   
    printf("Compiled with:");fprintf(ficlog,"Compiled with:");     printf("Compiled with:");fprintf(ficlog,"Compiled with:");
 #if defined(__clang__)  #if defined(__clang__)
    printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */     printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */
Line 5574  void syscompilerinfo() Line 5734  void syscompilerinfo()
    printf(" wtf-bit"); fprintf(ficlog," wtf-bit");/* wtf */     printf(" wtf-bit"); fprintf(ficlog," wtf-bit");/* wtf */
 #endif  #endif
   
 /* struct utsname sysInfo;  
   
    if (uname(&sysInfo) != -1) {  
      printf(" %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);  
      fprintf(ficlog," %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);  
    }  
    else  
       perror("uname() error");  
           */  
 #if defined(__GNUC__)  #if defined(__GNUC__)
 # if defined(__GNUC_PATCHLEVEL__)  # if defined(__GNUC_PATCHLEVEL__)
 #  define __GNUC_VERSION__ (__GNUC__ * 10000 \  #  define __GNUC_VERSION__ (__GNUC__ * 10000 \
Line 5594  void syscompilerinfo() Line 5745  void syscompilerinfo()
 # endif  # endif
    printf(" using GNU C version %d.\n", __GNUC_VERSION__);     printf(" using GNU C version %d.\n", __GNUC_VERSION__);
    fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__);     fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__);
    printf("GNU libc version: %s\n", gnu_get_libc_version());   
    fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version());   
   
      if (uname(&sysInfo) != -1) {
        printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
        fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
      }
      else
         perror("uname() error");
      //#ifndef __INTEL_COMPILER 
   #if !defined (__INTEL_COMPILER) && !defined(__APPLE__)
      printf("GNU libc version: %s\n", gnu_get_libc_version()); 
      fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version());
   #endif
 #endif  #endif
   
    //   void main()     //   void main()
    //   {     //   {
 #if defined(_MSC_VER)  #if defined(_MSC_VER)
    if (IsWow64()){     if (IsWow64()){
            printf("The program (probably compile for 32bit) is running under WOW64 (64bit) emulation.\n");             printf("The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
            fprintf(ficlog, "The program (ie 32bit) is running under WOW64 (64bit) emulation.\n");             fprintf(ficlog, "The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
    }     }
    else{     else{
            printf("The process is not running under WOW64 (i.e probably on a 64bits windows).\n");             printf("The process is not running under WOW64 (i.e probably on a 64bit Windows).\n");
            frintf(ficlog,"The programm is not running under WOW64 (i.e probably on a 64bits windows).\n");             fprintf(ficlog,"The programm is not running under WOW64 (i.e probably on a 64bit Windows).\n");
    }     }
    //      printf("\nPress Enter to continue...");     //      printf("\nPress Enter to continue...");
    //      getchar();     //      getchar();
Line 5619  void syscompilerinfo() Line 5779  void syscompilerinfo()
   
  }   }
   
   int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
     int i, j, k, i1 ;
     double ftolpl = 1.e-10;
     double age, agebase, agelim;
   
       strcpy(filerespl,"pl");
       strcat(filerespl,fileres);
       if((ficrespl=fopen(filerespl,"w"))==NULL) {
         printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
         fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
       }
       printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
       fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
       pstamp(ficrespl);
       fprintf(ficrespl,"# Period (stable) prevalence \n");
       fprintf(ficrespl,"#Age ");
       for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
       fprintf(ficrespl,"\n");
     
       /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
   
       agebase=ageminpar;
       agelim=agemaxpar;
   
       i1=pow(2,cptcoveff);
       if (cptcovn < 1){i1=1;}
   
       for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
         //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
           k=k+1;
           /* to clean */
           //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]);
           fprintf(ficrespl,"\n#******");
           printf("\n#******");
           fprintf(ficlog,"\n#******");
           for(j=1;j<=cptcoveff;j++) {
             fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           }
           fprintf(ficrespl,"******\n");
           printf("******\n");
           fprintf(ficlog,"******\n");
   
           fprintf(ficrespl,"#Age ");
           for(j=1;j<=cptcoveff;j++) {
             fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           }
           for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
           fprintf(ficrespl,"\n");
           
           for (age=agebase; age<=agelim; age++){
           /* for (age=agebase; age<=agebase; age++){ */
             prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
             fprintf(ficrespl,"%.0f ",age );
             for(j=1;j<=cptcoveff;j++)
               fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             for(i=1; i<=nlstate;i++)
               fprintf(ficrespl," %.5f", prlim[i][i]);
             fprintf(ficrespl,"\n");
           } /* Age */
           /* was end of cptcod */
       } /* cptcov */
           return 0;
   }
   
   int hPijx(double *p, int bage, int fage){
       /*------------- h Pij x at various ages ------------*/
   
     int stepsize;
     int agelim;
     int hstepm;
     int nhstepm;
     int h, i, i1, j, k;
   
     double agedeb;
     double ***p3mat;
   
       strcpy(filerespij,"pij");  strcat(filerespij,fileres);
       if((ficrespij=fopen(filerespij,"w"))==NULL) {
         printf("Problem with Pij resultfile: %s\n", filerespij); return 1;
         fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1;
       }
       printf("Computing pij: result on file '%s' \n", filerespij);
       fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);
     
       stepsize=(int) (stepm+YEARM-1)/YEARM;
       /*if (stepm<=24) stepsize=2;*/
   
       agelim=AGESUP;
       hstepm=stepsize*YEARM; /* Every year of age */
       hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
   
       /* hstepm=1;   aff par mois*/
       pstamp(ficrespij);
       fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
       i1= pow(2,cptcoveff);
      /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
      /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
      /*   k=k+1;  */
       for (k=1; k <= (int) pow(2,cptcoveff); k++){
         fprintf(ficrespij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++) 
           fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficrespij,"******\n");
         
         for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
           nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
           
           /*        nhstepm=nhstepm*YEARM; aff par mois*/
           
           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           oldm=oldms;savm=savms;
           hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
           fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
           for(i=1; i<=nlstate;i++)
             for(j=1; j<=nlstate+ndeath;j++)
               fprintf(ficrespij," %1d-%1d",i,j);
           fprintf(ficrespij,"\n");
           for (h=0; h<=nhstepm; h++){
             /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
             fprintf(ficrespij,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm );
             for(i=1; i<=nlstate;i++)
               for(j=1; j<=nlstate+ndeath;j++)
                 fprintf(ficrespij," %.5f", p3mat[i][j][h]);
             fprintf(ficrespij,"\n");
           }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           fprintf(ficrespij,"\n");
         }
         /*}*/
       }
           return 0;
   }
   
   
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
 /***********************************************/  /***********************************************/
Line 5726  int main(int argc, char *argv[]) Line 6025  int main(int argc, char *argv[])
   
   nberr=0; /* Number of errors and warnings */    nberr=0; /* Number of errors and warnings */
   nbwarn=0;    nbwarn=0;
   #ifdef WIN32
     _getcwd(pathcd, size);
   #else
   getcwd(pathcd, size);    getcwd(pathcd, size);
   #endif
   
   printf("\n%s\n%s",version,fullversion);    printf("\n%s\n%s",version,fullversion);
   if(argc <=1){    if(argc <=1){
Line 5762  int main(int argc, char *argv[]) Line 6065  int main(int argc, char *argv[])
   /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */    /* Split argv[1]=pathtot, parameter file name to get path, optionfile, extension and name */
   split(pathtot,path,optionfile,optionfilext,optionfilefiname);    split(pathtot,path,optionfile,optionfilext,optionfilefiname);
   printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);    printf("\npathtot=%s,\npath=%s,\noptionfile=%s \noptionfilext=%s \noptionfilefiname=%s\n",pathtot,path,optionfile,optionfilext,optionfilefiname);
   #ifdef WIN32
     _chdir(path); /* Can be a relative path */
     if(_getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */
   #else
   chdir(path); /* Can be a relative path */    chdir(path); /* Can be a relative path */
   if(getcwd(pathcd,MAXLINE) > 0) /* So pathcd is the full path */    if (getcwd(pathcd, MAXLINE) > 0) /* So pathcd is the full path */
     printf("Current directory %s!\n",pathcd);  #endif
     printf("Current directory %s!\n",pathcd);
   strcpy(command,"mkdir ");    strcpy(command,"mkdir ");
   strcat(command,optionfilefiname);    strcat(command,optionfilefiname);
   if((outcmd=system(command)) != 0){    if((outcmd=system(command)) != 0){
Line 6264  Title=%s <br>Datafile=%s Firstpass=%d La Line 6572  Title=%s <br>Datafile=%s Firstpass=%d La
   
   strcpy(pathr,path);    strcpy(pathr,path);
   strcat(pathr,optionfilefiname);    strcat(pathr,optionfilefiname);
   #ifdef WIN32
     _chdir(optionfilefiname); /* Move to directory named optionfile */
   #else
   chdir(optionfilefiname); /* Move to directory named optionfile */    chdir(optionfilefiname); /* Move to directory named optionfile */
   #endif
             
       
   /* Calculates basic frequencies. Computes observed prevalence at single age    /* Calculates basic frequencies. Computes observed prevalence at single age
      and prints on file fileres'p'. */       and prints on file fileres'p'. */
Line 6748  Interval (in months) between two waves: Line 7061  Interval (in months) between two waves:
   
   
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/      /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
 #include "prevlim.h"  /* Use ficrespl, ficlog */      /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */
       prlim=matrix(1,nlstate,1,nlstate);
       prevalence_limit(p, prlim,  ageminpar, agemaxpar);
     fclose(ficrespl);      fclose(ficrespl);
   
 #ifdef FREEEXIT2  #ifdef FREEEXIT2
Line 6756  Interval (in months) between two waves: Line 7071  Interval (in months) between two waves:
 #endif  #endif
   
     /*------------- h Pij x at various ages ------------*/      /*------------- h Pij x at various ages ------------*/
 #include "hpijx.h"      /*#include "hpijx.h"*/
       hPijx(p, bage, fage);
     fclose(ficrespij);      fclose(ficrespij);
   
   /*-------------- Variance of one-step probabilities---*/    /*-------------- Variance of one-step probabilities---*/
Line 7058  Interval (in months) between two waves: Line 7374  Interval (in months) between two waves:
   
   
    printf("Before Current directory %s!\n",pathcd);     printf("Before Current directory %s!\n",pathcd);
   #ifdef WIN32
      if (_chdir(pathcd) != 0)
              printf("Can't move to directory %s!\n",path);
      if(_getcwd(pathcd,MAXLINE) > 0)
   #else
    if(chdir(pathcd) != 0)     if(chdir(pathcd) != 0)
     printf("Can't move to directory %s!\n",path);             printf("Can't move to directory %s!\n", path);
   if(getcwd(pathcd,MAXLINE) > 0)     if (getcwd(pathcd, MAXLINE) > 0)
   #endif 
     printf("Current directory %s!\n",pathcd);      printf("Current directory %s!\n",pathcd);
   /*strcat(plotcmd,CHARSEPARATOR);*/    /*strcat(plotcmd,CHARSEPARATOR);*/
   sprintf(plotcmd,"gnuplot");    sprintf(plotcmd,"gnuplot");
Removed from v.1.175  
changed lines
  Added in v.1.184

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>