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

Diff for /imach/src/imach.c between versions 1.168 and 1.185

version 1.168, 2014/12/22 15:17:42 version 1.185, 2015/03/11 13:26:42
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.185  2015/03/11 13:26:42  brouard
     Summary: Inclusion of compile and links command line for Intel Compiler
   
     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
     *** empty log message ***
   
     Revision 1.174  2015/01/03 16:15:49  brouard
     Summary: Still in cross-compilation
   
     Revision 1.173  2015/01/03 12:06:26  brouard
     Summary: trying to detect cross-compilation
   
     Revision 1.172  2014/12/27 12:07:47  brouard
     Summary: Back from Visual Studio and Intel, options for compiling for Windows XP
   
     Revision 1.171  2014/12/23 13:26:59  brouard
     Summary: Back from Visual C
   
     Still problem with utsname.h on Windows
   
     Revision 1.170  2014/12/23 11:17:12  brouard
     Summary: Cleaning some \%% back to %%
   
     The escape was mandatory for a specific compiler (which one?), but too many warnings.
   
     Revision 1.169  2014/12/22 23:08:31  brouard
     Summary: 0.98p
   
     Outputs some informations on compiler used, OS etc. Testing on different platforms.
   
   Revision 1.168  2014/12/22 15:17:42  brouard    Revision 1.168  2014/12/22 15:17:42  brouard
   Summary: udate    Summary: update
   
   Revision 1.167  2014/12/22 13:50:56  brouard    Revision 1.167  2014/12/22 13:50:56  brouard
   Summary: Testing uname and compiler version and if compiled 32 or 64    Summary: Testing uname and compiler version and if compiled 32 or 64
Line 514 Line 579
 */  */
   
 #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 522 Line 589
   
 #ifdef _WIN32  #ifdef _WIN32
 #include <io.h>  #include <io.h>
   #include <windows.h>
   #include <tchar.h>
 #else  #else
 #include <unistd.h>  #include <unistd.h>
 #endif  #endif
   
 #include <limits.h>  #include <limits.h>
 #include <sys/types.h>  #include <sys/types.h>
 #include <sys/utsname.h>  
   #if defined(__GNUC__)
   #include <sys/utsname.h> /* Doesn't work on Windows */
   #endif
   
 #include <sys/stat.h>  #include <sys/stat.h>
 #include <errno.h>  #include <errno.h>
 /* extern int errno; */  /* extern int errno; */
Line 593  typedef struct { Line 666  typedef struct {
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.98nY, 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 709  static double maxarg1,maxarg2; Line 782  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 768  static int split( char *path, char *dirc Line 841  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 1225  double brent(double ax, double bx, doubl Line 1302  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 1248  double brent(double ax, double bx, doubl Line 1325  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 1256  void mnbrak(double *ax, double *bx, doub Line 1337  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 1275  void mnbrak(double *ax, double *bx, doub Line 1360  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 1356  void powell(double p[], double **xi, int Line 1493  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 1415  void powell(double p[], double **xi, int Line 1553  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 1436  void powell(double p[], double **xi, int Line 1579  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 1468  void powell(double p[], double **xi, int Line 1611  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 1498  void powell(double p[], double **xi, int Line 1642  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 :\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 1528  double **prevalim(double **prlim, int nl Line 1682  double **prevalim(double **prlim, int nl
 {  {
   /* Computes the prevalence limit in each live state at age x by left multiplying the unit    /* Computes the prevalence limit in each live state at age x by left multiplying the unit
      matrix by transitions matrix until convergence is reached */       matrix by transitions matrix until convergence is reached */
     
   int i, ii,j,k;    int i, ii,j,k;
   double min, max, maxmin, maxmax,sumnew=0.;    double min, max, maxmin, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **pmij();    double **out, cov[NCOVMAX+1], **pmij();
   double **newm;    double **newm;
   double agefin, delaymax=50 ; /* Max number of years to converge */    double agefin, delaymax=50 ; /* Max number of years to converge */
     
   for (ii=1;ii<=nlstate+ndeath;ii++)    for (ii=1;ii<=nlstate+ndeath;ii++)
     for (j=1;j<=nlstate+ndeath;j++){      for (j=1;j<=nlstate+ndeath;j++){
       oldm[ii][j]=(ii==j ? 1.0 : 0.0);        oldm[ii][j]=(ii==j ? 1.0 : 0.0);
     }      }
     
    cov[1]=1.;    cov[1]=1.;
      
  /* Even if hstepm = 1, at least one multiplication by the unit matrix */    /* Even if hstepm = 1, at least one multiplication by the unit matrix */
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){    for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
     newm=savm;      newm=savm;
     /* Covariates have to be included here again */      /* Covariates have to be included here again */
Line 1580  double **prevalim(double **prlim, int nl Line 1734  double **prevalim(double **prlim, int nl
       }        }
       maxmin=max-min;        maxmin=max-min;
       maxmax=FMAX(maxmax,maxmin);        maxmax=FMAX(maxmax,maxmin);
     }      } /* j loop */
     if(maxmax < ftolpl){      if(maxmax < ftolpl){
       return prlim;        return prlim;
     }      }
   }    } /* age loop */
     return prlim; /* should not reach here */
 }  }
   
 /*************** transition probabilities ***************/   /*************** transition probabilities ***************/ 
Line 1867  double func( double *x) Line 2022  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 1899  double func( double *x) Line 2073  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 2234  void mlikeli(FILE *ficres,double p[], in Line 2412  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 2533  void  freqsummary(char fileres[], int ia Line 2711  void  freqsummary(char fileres[], int ia
   
   first=1;    first=1;
   
   /* for(k1=1; k1<=j ; k1++){   /* Loop on covariates */    /* for(k1=1; k1<=j ; k1++){ */  /* Loop on covariates */
   /*  for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */    /*  for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */
   /*    j1++;    /*    j1++; */
 */  
   for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){    for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){
       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
         scanf("%d", i);*/          scanf("%d", i);*/
Line 2913  void tricode(int *Tvar, int **nbcode, in Line 3090  void tricode(int *Tvar, int **nbcode, in
 {  {
   /**< Uses cptcovn+2*cptcovprod as the number of covariates */    /**< Uses cptcovn+2*cptcovprod as the number of covariates */
   /*      Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1     /*      Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 
   /* Boring subroutine which should only output nbcode[Tvar[j]][k]     * Boring subroutine which should only output nbcode[Tvar[j]][k]
    * Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2)     * Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2)
   /* nbcode[Tvar[j]][1]=      * nbcode[Tvar[j]][1]= 
   */    */
   
   int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;    int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
Line 3343  void varevsij(char optionfilefiname[], d Line 3520  void varevsij(char optionfilefiname[], d
   /* Variance of health expectancies */    /* Variance of health expectancies */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
   /* double **newm;*/    /* double **newm;*/
     /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/
     
     int movingaverage();
   double **dnewm,**doldm;    double **dnewm,**doldm;
   double **dnewmp,**doldmp;    double **dnewmp,**doldmp;
   int i, j, nhstepm, hstepm, h, nstepm ;    int i, j, nhstepm, hstepm, h, nstepm ;
Line 3620  void varevsij(char optionfilefiname[], d Line 3800  void varevsij(char optionfilefiname[], d
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */  /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */  /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
   fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));    fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95\%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));    fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
   fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);    fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
Line 4041  To be simple, these graphs help to under Line 4221  To be simple, these graphs help to under
           } /* k12 */            } /* k12 */
         } /*l1 */          } /*l1 */
       }/* k1 */        }/* k1 */
       /* } /* loop covariates */        /* } */ /* loop covariates */
   }    }
   free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);    free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
   free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);    free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
Line 4218  void printinggnuplot(char fileres[], cha Line 4398  void printinggnuplot(char fileres[], cha
      fprintf(ficgp,"set xlabel \"Age\" \n\       fprintf(ficgp,"set xlabel \"Age\" \n\
 set ylabel \"Probability\" \n\  set ylabel \"Probability\" \n\
 set ter png small size 320, 240\n\  set ter png small size 320, 240\n\
 plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
   
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else        fprintf(ficgp," \%%*lf (\%%*lf)");         else        fprintf(ficgp," %%*lf (%%*lf)");
      }       }
      fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);       fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," %%*lf (%%*lf)");
      }        } 
      fprintf(ficgp,"\" t\"95\%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);        fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); 
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," %%*lf (%%*lf)");
      }         }  
      fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));       fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
    }     }
Line 4245  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 4425  plot [%.f:%.f] \"%s\" every :::%d::%d u
           
     for (i=1; i<= nlstate+1 ; i ++) {      for (i=1; i<= nlstate+1 ; i ++) {
       k=2*i;        k=2*i;
       fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);        fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
       for (j=1; j<= nlstate+1 ; j ++) {        for (j=1; j<= nlstate+1 ; j ++) {
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");          if (j==i) fprintf(ficgp," %%lf (%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");          else fprintf(ficgp," %%*lf (%%*lf)");
       }           }   
       if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");        if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
       else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);        else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
       fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);        fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
       for (j=1; j<= nlstate+1 ; j ++) {        for (j=1; j<= nlstate+1 ; j ++) {
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");          if (j==i) fprintf(ficgp," %%lf (%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");          else fprintf(ficgp," %%*lf (%%*lf)");
       }           }   
       fprintf(ficgp,"\" t\"\" w l lt 0,");        fprintf(ficgp,"\" t\"\" w l lt 0,");
       fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);        fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
       for (j=1; j<= nlstate+1 ; j ++) {        for (j=1; j<= nlstate+1 ; j ++) {
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");          if (j==i) fprintf(ficgp," %%lf (%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");          else fprintf(ficgp," %%*lf (%%*lf)");
       }           }   
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");        if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
       else fprintf(ficgp,"\" t\"\" w l lt 0,");        else fprintf(ficgp,"\" t\"\" w l lt 0,");
Line 4429  int movingaverage(double ***probs, doubl Line 4609  int movingaverage(double ***probs, doubl
   
   
 /************** Forecasting ******************/  /************** Forecasting ******************/
 prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){  void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
   /* proj1, year, month, day of starting projection     /* proj1, year, month, day of starting projection 
      agemin, agemax range of age       agemin, agemax range of age
      dateprev1 dateprev2 range of dates during which prevalence is computed       dateprev1 dateprev2 range of dates during which prevalence is computed
Line 4552  prevforecast(char fileres[], double anpr Line 4732  prevforecast(char fileres[], double anpr
 }  }
   
 /************** Forecasting *****not tested NB*************/  /************** Forecasting *****not tested NB*************/
 populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){  void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
       
   int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;    int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
   int *popage;    int *popage;
Line 5028  int readdata(char datafile[], int firsto Line 5208  int readdata(char datafile[], int firsto
               
       strcpy(line,stra);        strcpy(line,stra);
       cutv(stra, strb,line,' ');        cutv(stra, strb,line,' ');
       if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){        if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
       }        }
       else  if(iout=sscanf(strb,"%s.",dummy) != 0){        else  if( (iout=sscanf(strb,"%s.",dummy)) != 0){
         month=99;          month=99;
         year=9999;          year=9999;
       }else{        }else{
Line 5044  int readdata(char datafile[], int firsto Line 5224  int readdata(char datafile[], int firsto
     } /* ENd Waves */      } /* ENd Waves */
           
     cutv(stra, strb,line,' ');       cutv(stra, strb,line,' '); 
     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){      if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
     }      }
     else  if(iout=sscanf(strb,"%s.",dummy) != 0){      else  if( (iout=sscanf(strb,"%s.",dummy)) != 0){
       month=99;        month=99;
       year=9999;        year=9999;
     }else{      }else{
Line 5059  int readdata(char datafile[], int firsto Line 5239  int readdata(char datafile[], int firsto
     strcpy(line,stra);      strcpy(line,stra);
           
     cutv(stra, strb,line,' ');       cutv(stra, strb,line,' '); 
     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){      if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
     }      }
     else  if(iout=sscanf(strb,"%s.", dummy) != 0){      else  if( (iout=sscanf(strb,"%s.", dummy)) != 0){
       month=99;        month=99;
       year=9999;        year=9999;
     }else{      }else{
Line 5159  void removespace(char *str) { Line 5339  void removespace(char *str) {
   do    do
     while (*p2 == ' ')      while (*p2 == ' ')
       p2++;        p2++;
   while (*p1++ = *p2++);    while (*p1++ == *p2++);
 }  }
   
 int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:  int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:
Line 5189  int decodemodel ( char model[], int last Line 5369  int decodemodel ( char model[], int last
     cptcovs=j+1-j1; /**<  Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2  */      cptcovs=j+1-j1; /**<  Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2  */
     cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/      cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/
                   /* including age products which are counted in cptcovage.                    /* including age products which are counted in cptcovage.
                   /* but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */                    * but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */
     cptcovprod=j1; /**< Number of products  V1*V2 +v3*age = 2 */      cptcovprod=j1; /**< Number of products  V1*V2 +v3*age = 2 */
     cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1  */      cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1  */
     strcpy(modelsav,model);       strcpy(modelsav,model); 
Line 5331  int decodemodel ( char model[], int last Line 5511  int decodemodel ( char model[], int last
     return (1);      return (1);
 }  }
   
 calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )  int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
 {  {
   int i, m;    int i, m;
   
Line 5342  calandcheckages(int imx, int maxwav, dou Line 5522  calandcheckages(int imx, int maxwav, dou
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){        if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
         *nberr++;          *nberr = *nberr + 1;
         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);          printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);          fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){        if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
         *nberr++;          (*nberr)++;
         printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]);           printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
         fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]);           fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
         s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */          s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
Line 5361  calandcheckages(int imx, int maxwav, dou Line 5541  calandcheckages(int imx, int maxwav, dou
     for(m=firstpass; (m<= lastpass); m++){      for(m=firstpass; (m<= lastpass); m++){
       if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){        if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
         if (s[m][i] >= nlstate+1) {          if (s[m][i] >= nlstate+1) {
           if(agedc[i]>0)            if(agedc[i]>0){
             if((int)moisdc[i]!=99 && (int)andc[i]!=9999)              if((int)moisdc[i]!=99 && (int)andc[i]!=9999){
               agev[m][i]=agedc[i];                agev[m][i]=agedc[i];
           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/            /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
             else {              }else {
               if ((int)andc[i]!=9999){                if ((int)andc[i]!=9999){
                 nbwarn++;                  nbwarn++;
                 printf("Warning negative age at death: %ld line:%d\n",num[i],i);                  printf("Warning negative age at death: %ld line:%d\n",num[i],i);
Line 5373  calandcheckages(int imx, int maxwav, dou Line 5553  calandcheckages(int imx, int maxwav, dou
                 agev[m][i]=-1;                  agev[m][i]=-1;
               }                }
             }              }
             } /* agedc > 0 */
         }          }
         else if(s[m][i] !=9){ /* Standard case, age in fractional          else if(s[m][i] !=9){ /* Standard case, age in fractional
                                  years but with the precision of a month */                                   years but with the precision of a month */
Line 5403  calandcheckages(int imx, int maxwav, dou Line 5584  calandcheckages(int imx, int maxwav, dou
   for (i=1; i<=imx; i++)  {    for (i=1; i<=imx; i++)  {
     for(m=firstpass; (m<=lastpass); m++){      for(m=firstpass; (m<=lastpass); m++){
       if (s[m][i] > (nlstate+ndeath)) {        if (s[m][i] > (nlstate+ndeath)) {
         *nberr++;          (*nberr)++;
         printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);               printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
         fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);               fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
         return 1;          return 1;
Line 5428  calandcheckages(int imx, int maxwav, dou Line 5609  calandcheckages(int imx, int maxwav, dou
     return (1);      return (1);
 }  }
   
 syscompilerinfo()  #if defined(_MSC_VER)
   /*printf("Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/
   /*fprintf(ficlog, "Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/
   //#include "stdafx.h"
   //#include <stdio.h>
   //#include <tchar.h>
   //#include <windows.h>
   //#include <iostream>
   typedef BOOL(WINAPI *LPFN_ISWOW64PROCESS) (HANDLE, PBOOL);
   
   LPFN_ISWOW64PROCESS fnIsWow64Process;
   
   BOOL IsWow64()
   {
           BOOL bIsWow64 = FALSE;
   
           //typedef BOOL (APIENTRY *LPFN_ISWOW64PROCESS)
           //  (HANDLE, PBOOL);
   
           //LPFN_ISWOW64PROCESS fnIsWow64Process;
   
           HMODULE module = GetModuleHandle(_T("kernel32"));
           const char funcName[] = "IsWow64Process";
           fnIsWow64Process = (LPFN_ISWOW64PROCESS)
                   GetProcAddress(module, funcName);
   
           if (NULL != fnIsWow64Process)
           {
                   if (!fnIsWow64Process(GetCurrentProcess(),
                           &bIsWow64))
                           //throw std::exception("Unknown error");
                           printf("Unknown error\n");
           }
           return bIsWow64 != FALSE;
   }
   #endif
   
   void syscompilerinfo()
  {   {
    /* #include "syscompilerinfo.h"*/     /* #include "syscompilerinfo.h"*/
 #include <gnu/libc-version.h>     /* command line Intel compiler 32bit windows, XP compatible:*/
      /* /GS /W3 /Gy
         /Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D
         "_CONSOLE" /D "_LIB" /D "_USING_V110_SDK71_" /D "_UNICODE" /D
         "UNICODE" /Qipo /Zc:forScope /Gd /Oi /MT /Fa"Release\" /EHsc /nologo
         /Fo"Release\" /Qprof-dir "Release\" /Fp"Release\IMaCh.pch" 
         /* 64 bits */
      /*
        /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" */
      /* Optimization are useless and O3 is slower than O2 */
      /*
        /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" 
      */
      /* Link is $/ /* /OUT:"visual studio
         2013\Projects\IMaCh\Release\IMaCh.exe" /MANIFEST /NXCOMPAT
         /PDB:"visual studio
         2013\Projects\IMaCh\Release\IMaCh.pdb" /DYNAMICBASE
         "kernel32.lib" "user32.lib" "gdi32.lib" "winspool.lib"
         "comdlg32.lib" "advapi32.lib" "shell32.lib" "ole32.lib"
         "oleaut32.lib" "uuid.lib" "odbc32.lib" "odbccp32.lib"
         /MACHINE:X86 /OPT:REF /SAFESEH /INCREMENTAL:NO
         /SUBSYSTEM:CONSOLE",5.01" /MANIFESTUAC:"level='asInvoker'
         uiAccess='false'"
         /ManifestFile:"Release\IMaCh.exe.intermediate.manifest" /OPT:ICF
         /NOLOGO /TLBID:1
      */
   #if defined __INTEL_COMPILER
 #if defined(__GNUC__)  #if defined(__GNUC__)
 # if defined(__GNUC_PATCHLEVEL__)          struct utsname sysInfo;  /* For Intel on Linux and OS/X */
 #  define __GNUC_VERSION__ (__GNUC__ * 10000 \  #endif
                             + __GNUC_MINOR__ * 100 \  #elif defined(__GNUC__) 
                             + __GNUC_PATCHLEVEL__)  #ifndef  __APPLE__
 # else  #include <gnu/libc-version.h>  /* Only on gnu */
 #  define __GNUC_VERSION__ (__GNUC__ * 10000 \  
                             + __GNUC_MINOR__ * 100)  
 # endif  
 #endif  #endif
      struct utsname sysInfo;
      int cross = CROSS;
      if (cross){
              printf("Cross-");
              fprintf(ficlog, "Cross-");
      }
   #endif
   
   #include <stdint.h>
   
      printf("Compiled with:");fprintf(ficlog,"Compiled with:");
   #if defined(__clang__)
      printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */
   #endif
   #if defined(__ICC) || defined(__INTEL_COMPILER)
      printf(" Intel ICC/ICPC");fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */
   #endif
   #if defined(__GNUC__) || defined(__GNUG__)
      printf(" GNU GCC/G++");fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */
   #endif
   #if defined(__HP_cc) || defined(__HP_aCC)
      printf(" Hewlett-Packard C/aC++");fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */
   #endif
   #if defined(__IBMC__) || defined(__IBMCPP__)
      printf(" IBM XL C/C++"); fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */
   #endif
   #if defined(_MSC_VER)
      printf(" Microsoft Visual Studio");fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */
   #endif
   #if defined(__PGI)
      printf(" Portland Group PGCC/PGCPP");fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */
   #endif
   #if defined(__SUNPRO_C) || defined(__SUNPRO_CC)
      printf(" Oracle Solaris Studio");fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */
   #endif
      printf(" for ");fprintf(ficlog," for ");
      
 // http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros  // http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros
 #ifdef _WIN32 // note the underscore: without it, it's not msdn official!  #ifdef _WIN32 // note the underscore: without it, it's not msdn official!
     // Windows (x64 and x86)      // Windows (x64 and x86)
      printf("Windows (x64 and x86) ");fprintf(ficlog,"Windows (x64 and x86) ");
 #elif __unix__ // all unices, not all compilers  #elif __unix__ // all unices, not all compilers
     // Unix      // Unix
      printf("Unix ");fprintf(ficlog,"Unix ");
 #elif __linux__  #elif __linux__
     // linux      // linux
      printf("linux ");fprintf(ficlog,"linux ");
 #elif __APPLE__  #elif __APPLE__
     // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though...      // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though..
      printf("Mac OS ");fprintf(ficlog,"Mac OS ");
 #endif  #endif
   
 /*  __MINGW32__   */  /*  __MINGW32__   */
Line 5463  syscompilerinfo() Line 5751  syscompilerinfo()
 /* _WIN64  // Defined for applications for Win64. */  /* _WIN64  // Defined for applications for Win64. */
 /* _M_X64 // Defined for compilations that target x64 processors. */  /* _M_X64 // Defined for compilations that target x64 processors. */
 /* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */  /* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */
 #include <stdint.h>  
 #if UINTPTR_MAX == 0xffffffff  #if UINTPTR_MAX == 0xffffffff
    printf("32-bit \n"); /* 32-bit */     printf(" 32-bit"); fprintf(ficlog," 32-bit");/* 32-bit */
 #elif UINTPTR_MAX == 0xffffffffffffffff  #elif UINTPTR_MAX == 0xffffffffffffffff
   printf("64-bit \n");/* 64-bit */     printf(" 64-bit"); fprintf(ficlog," 64-bit");/* 64-bit */
 #else  #else
  printf("wtf-bit \n"); /* wtf */     printf(" wtf-bit"); fprintf(ficlog," wtf-bit");/* wtf */
 #endif  #endif
   
 struct utsname sysInfo;  #if defined(__GNUC__)
   # if defined(__GNUC_PATCHLEVEL__)
   #  define __GNUC_VERSION__ (__GNUC__ * 10000 \
                               + __GNUC_MINOR__ * 100 \
                               + __GNUC_PATCHLEVEL__)
   # else
   #  define __GNUC_VERSION__ (__GNUC__ * 10000 \
                               + __GNUC_MINOR__ * 100)
   # endif
      printf(" using GNU C version %d.\n", __GNUC_VERSION__);
      fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__);
   
    if (uname(&sysInfo) != -1) {     if (uname(&sysInfo) != -1) {
       puts(sysInfo.sysname);       printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
       puts(sysInfo.nodename);       fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
       puts(sysInfo.release);  
       puts(sysInfo.version);  
       puts(sysInfo.machine);  
    }     }
    else     else
       perror("uname() error");        perror("uname() error");
    printf("GNU C version %d\n", __GNUC_VERSION__);     //#ifndef __INTEL_COMPILER 
   printf("GNU libc version: %s\n", gnu_get_libc_version());  #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
   
      //   void main()
      //   {
   #if defined(_MSC_VER)
      if (IsWow64()){
              printf("The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
              fprintf(ficlog, "The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
      }
      else{
              printf("The process is not running under WOW64 (i.e probably on a 64bit Windows).\n");
              fprintf(ficlog,"The programm is not running under WOW64 (i.e probably on a 64bit Windows).\n");
      }
      //      printf("\nPress Enter to continue...");
      //      getchar();
      //   }
   
   #endif
      
   
  }   }
   
   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 5595  int main(int argc, char *argv[]) Line 6051  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 5631  int main(int argc, char *argv[]) Line 6091  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){
     printf("Problem creating directory or it already exists %s%s, err=%d\n",path,optionfilefiname,outcmd);      printf("Directory already exists (or can't create it) %s%s, err=%d\n",path,optionfilefiname,outcmd);
     /* fprintf(ficlog,"Problem creating directory %s%s\n",path,optionfilefiname); */      /* fprintf(ficlog,"Problem creating directory %s%s\n",path,optionfilefiname); */
     /* fclose(ficlog); */      /* fclose(ficlog); */
 /*     exit(1); */  /*     exit(1); */
Line 6133  Title=%s <br>Datafile=%s Firstpass=%d La Line 6598  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 6617  Interval (in months) between two waves: Line 7087  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 6625  Interval (in months) between two waves: Line 7097  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 6927  Interval (in months) between two waves: Line 7400  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.168  
changed lines
  Added in v.1.185

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