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

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

version 1.182, 2015/02/12 08:19:57 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    Revision 1.182  2015/02/12 08:19:57  brouard
   Summary: Trying to keep directest which seems simpler and more general    Summary: Trying to keep directest which seems simpler and more general
   Author: Nicolas Brouard    Author: Nicolas Brouard
Line 565 Line 576
 */  */
   
 #define POWELL /* Instead of NLOPT */  #define POWELL /* Instead of NLOPT */
 #define POWELLDIRECT /* Directest to decide new direction instead of Powell test */  /* #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 651  typedef struct { Line 663  typedef struct {
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.98p, Février 2015,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 767  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 826  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 1283  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 1306  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 1314  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 1333  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 1474  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) { /* We are keeping the max gain on each of the n directions         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                                         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.                                        in comparison with the 'probable' gain, mu^2, with the last average direction.
Line 1546  void powell(double p[], double **xi, int Line 1622  void powell(double p[], double **xi, int
       /* 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 */
       /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */        /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < 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); /* Intel compiler doesn't work on one line */        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);
   #endif
       directest = fp-2.0*(*fret)+fptt - 2.0 * del; /* If del was big enough we change it for a new direction */        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);        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);
Line 1560  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
 #ifdef POWELLDIRECT  #ifdef POWELLORIGINAL
         if (t < 0.0) { /* Then we use it for new direction */
   #else
       if (directest*t < 0.0) { /* Contradiction between both tests */        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);        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);        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);        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);        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 */        if (directest < 0.0) { /* Then we use it for new direction */
 #else  
       if (t < 0.0) { /* Then we use it for new direction */  
 #endif  #endif
         linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/          linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
Line 1940  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 1972  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 5544  BOOL IsWow64() Line 5646  BOOL IsWow64()
 void syscompilerinfo()  void syscompilerinfo()
  {   {
    /* #include "syscompilerinfo.h"*/     /* #include "syscompilerinfo.h"*/
           /* 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 __INTEL_COMPILER
 #if defined(__GNUC__)  #if defined(__GNUC__)
         struct utsname sysInfo;  /* For Intel on Linux and OS/X */          struct utsname sysInfo;  /* For Intel on Linux and OS/X */
Line 5733  int prevalence_limit(double *p, double * Line 5844  int prevalence_limit(double *p, double *
         } /* Age */          } /* Age */
         /* was end of cptcod */          /* was end of cptcod */
     } /* cptcov */      } /* cptcov */
           return 0;
 }  }
   
 int hPijx(double *p, int bage, int fage){  int hPijx(double *p, int bage, int fage){
Line 5766  int hPijx(double *p, int bage, int fage) Line 5878  int hPijx(double *p, int bage, int fage)
     pstamp(ficrespij);      pstamp(ficrespij);
     fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");      fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
     i1= pow(2,cptcoveff);      i1= pow(2,cptcoveff);
    for(cptcov=1,k=0;cptcov<=i1;cptcov++){     /* for(cptcov=1,k=0;cptcov<=i1;cptcov++){ */
       /*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/     /*    /\*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*\/ */
         k=k+1;      /*   k=k+1;  */
     /* for (k=1; k <= (int) pow(2,cptcoveff); k++){*/      for (k=1; k <= (int) pow(2,cptcoveff); k++){
         fprintf(ficrespij,"\n#****** ");        fprintf(ficrespij,"\n#****** ");
         for(j=1;j<=cptcoveff;j++)         for(j=1;j<=cptcoveff;j++) 
           fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);          fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
         fprintf(ficrespij,"******\n");        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 */
                   
         for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */          /*        nhstepm=nhstepm*YEARM; aff par mois*/
           nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */           
           nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */          p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           oldm=oldms;savm=savms;
           /*      nhstepm=nhstepm*YEARM; aff par mois*/          hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
           fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
           p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          for(i=1; i<=nlstate;i++)
           oldm=oldms;savm=savms;            for(j=1; j<=nlstate+ndeath;j++)
           hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);                fprintf(ficrespij," %1d-%1d",i,j);
           fprintf(ficrespij,"# Cov Agex agex+h hpijx with 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(i=1; i<=nlstate;i++)
             for(j=1; j<=nlstate+ndeath;j++)              for(j=1; j<=nlstate+ndeath;j++)
               fprintf(ficrespij," %1d-%1d",i,j);                fprintf(ficrespij," %.5f", p3mat[i][j][h]);
           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");            fprintf(ficrespij,"\n");
         }          }
           free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
           fprintf(ficrespij,"\n");
         }
       /*}*/        /*}*/
     }      }
           return 0;
 }  }
   
   
Line 5912  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 5948  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 6450  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 7247  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.182  
changed lines
  Added in v.1.184

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