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

Diff for /imach/src/imach-0.98r.c between versions 1.1 and 1.2

version 1.1, 2016/02/19 09:46:04 version 1.2, 2019/05/21 15:55:06
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.2  2019/05/21 15:55:06  brouard
     Summary: RPM and SRPMS
   
   Revision 1.1  2016/02/19 09:46:04  brouard    Revision 1.1  2016/02/19 09:46:04  brouard
   Summary: Kind of 0.98r? series, starting with r7    Summary: Kind of 0.98r? series, starting with r7
   
Line 747  Back prevalence and projections: Line 750  Back prevalence and projections:
 #define DEBUGHESSIJ  #define DEBUGHESSIJ
 /* #define LINMINORIGINAL  /\* Don't use loop on scale in linmin (accepting nan)*\/ */  /* #define LINMINORIGINAL  /\* Don't use loop on scale in linmin (accepting nan)*\/ */
 #define POWELL /* Instead of NLOPT */  #define POWELL /* Instead of NLOPT */
 #define POWELLF1F3 /* Skip test */  /* #ifndef POWELLNOF3INFF1TEST  */
   /* #define POWELLNOF3INFF1TEST /\* Skip test *\/ */
   /* #endif */
 /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */  /* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
 /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */  /* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
   
Line 840  typedef struct { Line 845  typedef struct {
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
 #include "version.h"  #include "version.h"
 char version[]=__IMACH_VERSION__;  char version[]=__IMACH_VERSION_98R__;
 char copyright[]="October 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";  char copyright[]="February 2016,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015-2018";
 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 1545  double brent(double ax, double bx, doubl Line 1550  double brent(double ax, double bx, doubl
       etemp=e;         etemp=e; 
       e=d;         e=d; 
       if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))         if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) 
         d=CGOLD*(e=(x >= xm ? a-x : b-x));                                   d=CGOLD*(e=(x >= xm ? a-x : b-x)); 
       else {         else { 
         d=p/q;                                   d=p/q; 
         u=x+d;                                   u=x+d; 
         if (u-a < tol2 || b-u < tol2)                                   if (u-a < tol2 || b-u < tol2) 
           d=SIGN(tol1,xm-x);                                           d=SIGN(tol1,xm-x); 
       }         } 
     } else {       } else { 
       d=CGOLD*(e=(x >= xm ? a-x : b-x));         d=CGOLD*(e=(x >= xm ? a-x : b-x)); 
Line 1564  double brent(double ax, double bx, doubl Line 1569  double brent(double ax, double bx, doubl
     } 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; 
       }         } 
     }       } 
   }     } 
Line 1611  values at the three points, fa, fb , and Line 1616  values at the three points, fa, fb , and
   *cx=(*bx)+GOLD*(*bx-*ax);     *cx=(*bx)+GOLD*(*bx-*ax); 
   *fc=(*func)(*cx);     *fc=(*func)(*cx); 
 #ifdef DEBUG  #ifdef DEBUG
   printf("mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);    printf("mnbrak0 a=%lf *fa=%lf, b=%lf *fb=%lf, c=%lf *fc=%lf\n",*ax,*fa,*bx,*fb,*cx, *fc);
   fprintf(ficlog,"mnbrak0 *fb=%.12e *fc=%.12e\n",*fb,*fc);    fprintf(ficlog,"mnbrak0 a=%lf *fa=%lf, b=%lf *fb=%lf, c=%lf *fc=%lf\n",*ax,*fa,*bx,*fb,*cx, *fc);
 #endif  #endif
   while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc */    while (*fb > *fc) { /* Declining a,b,c with fa> fb > fc. If fc=inf it exits*/
     r=(*bx-*ax)*(*fb-*fc);       r=(*bx-*ax)*(*fb-*fc); 
     q=(*bx-*cx)*(*fb-*fa);       q=(*bx-*cx)*(*fb-*fa); /* What if fa=inf */
     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 abscissa 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 abscissa where function should be evaluated */      ulim=(*bx)+GLIMIT*(*cx-*bx); /* Maximum abscissa where function should be evaluated */
     if ((*bx-u)*(u-*cx) > 0.0) { /* if u_p is 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) */
       double A, fparabu;         double A, fparabu; 
       A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u);        A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u);
       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("\nmnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f, q=%lf < %lf=r)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu,q,r);
       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,"\nmnbrak (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf),  (*u=%.12f, fu=%.12lf, fparabu=%.12f, q=%lf < %lf=r)\n",*ax,*fa,*bx,*fb,*cx,*fc,u,fu, fparabu,q,r);
       /* And thus,it can be that fu > *fc even if fparabu < *fc */        /* 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),        /* 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) */          (*cx=10.098840694817, *fc=298946.631474258087),  (*u=9.852501168332, fu=298948.773013752128, fparabu=298945.434711494134) */
Line 1661  values at the three points, fa, fb , and Line 1666  values at the three points, fa, fb , and
 /*      fu = *fc; */  /*      fu = *fc; */
 /*      *fc =dum; */  /*      *fc =dum; */
 /*       } */  /*       } */
 #ifdef DEBUG  #ifdef DEBUGMNBRAK
       printf("mnbrak34  fu < or >= fc \n");                   double A, fparabu; 
       fprintf(ficlog, "mnbrak34 fu < fc\n");       A= (*fb - *fa)/(*bx-*ax)/(*bx+*ax-2*u);
        fparabu= *fa - A*(*ax-u)*(*ax-u);
        printf("\nmnbrak35 ax=%lf fa=%lf bx=%lf fb=%lf, u=%lf fp=%lf fu=%lf < or >= fc=%lf cx=%lf, q=%lf < %lf=r \n",*ax, *fa, *bx,*fb,u,fparabu,fu,*fc,*cx,q,r);
        fprintf(ficlog,"\nmnbrak35 ax=%lf fa=%lf bx=%lf fb=%lf, u=%lf fp=%lf fu=%lf < or >= fc=%lf cx=%lf, q=%lf < %lf=r \n",*ax, *fa, *bx,*fb,u,fparabu,fu,*fc,*cx,q,r);
 #endif  #endif
       dum=u; /* Shifting c and u */        dum=u; /* Shifting c and u */
       u = *cx;        u = *cx;
Line 1674  values at the three points, fa, fb , and Line 1682  values at the three points, fa, fb , and
 #endif  #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  #ifdef DEBUG
       printf("mnbrak2  u after c but before ulim\n");        printf("\nmnbrak2  u=%lf after c=%lf but before ulim\n",u,*cx);
       fprintf(ficlog, "mnbrak2 u after c but before ulim\n");        fprintf(ficlog,"\nmnbrak2  u=%lf after c=%lf but before ulim\n",u,*cx);
 #endif  #endif
       fu=(*func)(u);         fu=(*func)(u); 
       if (fu < *fc) {         if (fu < *fc) { 
 #ifdef DEBUG  #ifdef DEBUG
       printf("mnbrak2  u after c but before ulim AND fu < fc\n");                                  printf("\nmnbrak2  u=%lf after c=%lf but before ulim=%lf AND fu=%lf < %lf=fc\n",u,*cx,ulim,fu, *fc);
       fprintf(ficlog, "mnbrak2 u after c but before ulim AND fu <fc \n");                            fprintf(ficlog,"\nmnbrak2  u=%lf after c=%lf but before ulim=%lf AND fu=%lf < %lf=fc\n",u,*cx,ulim,fu, *fc);
   #endif
                             SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) 
                                   SHFT(*fb,*fc,fu,(*func)(u)) 
   #ifdef DEBUG
                                           printf("\nmnbrak2 shift GOLD c=%lf",*cx+GOLD*(*cx-*bx));
 #endif  #endif
         SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))   
         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  #ifdef DEBUG
       printf("mnbrak2  u outside ulim (verifying that ulim is beyond c)\n");        printf("\nmnbrak2  u=%lf outside ulim=%lf (verifying that ulim is beyond c=%lf)\n",u,ulim,*cx);
       fprintf(ficlog, "mnbrak2 u outside ulim (verifying that ulim is beyond c)\n");        fprintf(ficlog,"\nmnbrak2  u=%lf outside ulim=%lf (verifying that ulim is beyond c=%lf)\n",u,ulim,*cx);
 #endif  #endif
       u=ulim;         u=ulim; 
       fu=(*func)(u);         fu=(*func)(u); 
     } else { /* u could be left to b (if r > q parabola has a maximum) */      } else { /* u could be left to b (if r > q parabola has a maximum) */
 #ifdef DEBUG  #ifdef DEBUG
       printf("mnbrak2  u could be left to b (if r > q parabola has a maximum)\n");        printf("\nmnbrak2  u=%lf could be left to b=%lf (if r=%lf > q=%lf parabola has a maximum)\n",u,*bx,r,q);
       fprintf(ficlog, "mnbrak2  u could be left to b (if r > q parabola has a maximum)\n");        fprintf(ficlog,"\nmnbrak2  u=%lf could be left to b=%lf (if r=%lf > q=%lf parabola has a maximum)\n",u,*bx,r,q);
 #endif  #endif
       u=(*cx)+GOLD*(*cx-*bx);         u=(*cx)+GOLD*(*cx-*bx); 
       fu=(*func)(u);         fu=(*func)(u); 
   #ifdef DEBUG
         printf("\nmnbrak2 new u=%lf fu=%lf shifted gold left from c=%lf and b=%lf \n",u,fu,*cx,*bx);
         fprintf(ficlog,"\nmnbrak2 new u=%lf fu=%lf shifted gold left from c=%lf and b=%lf \n",u,fu,*cx,*bx);
   #endif
     } /* end tests */      } /* end tests */
     SHFT(*ax,*bx,*cx,u)       SHFT(*ax,*bx,*cx,u) 
     SHFT(*fa,*fb,*fc,fu)       SHFT(*fa,*fb,*fc,fu) 
 #ifdef DEBUG  #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);        printf("\nmnbrak2 shift (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc);
       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);        fprintf(ficlog, "\nmnbrak2 shift (*ax=%.12f, *fa=%.12lf), (*bx=%.12f, *fb=%.12lf), (*cx=%.12f, *fc=%.12lf)\n",*ax,*fa,*bx,*fb,*cx,*fc);
 #endif  #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) */    } /* end while; ie return (a, b, c, fa, fb, fc) such that a < b < c with f(a) > f(b) and fb < f(c) */
 }   } 
Line 1720  int ncom; Line 1735  int ncom;
 double *pcom,*xicom;  double *pcom,*xicom;
 double (*nrfunc)(double []);   double (*nrfunc)(double []); 
     
   #ifdef LINMINORIGINAL
 void linmin(double p[], double xi[], int n, double *fret,double (*func)(double []))   void linmin(double p[], double xi[], int n, double *fret,double (*func)(double [])) 
   #else
   void linmin(double p[], double xi[], int n, double *fret,double (*func)(double []), int *flat) 
   #endif
 {   { 
   double brent(double ax, double bx, double cx,     double brent(double ax, double bx, double cx, 
                double (*f)(double), double tol, double *xmin);                  double (*f)(double), double tol, double *xmin); 
Line 1764  void linmin(double p[], double xi[], int Line 1783  void linmin(double p[], double xi[], int
 #ifdef LINMINORIGINAL  #ifdef LINMINORIGINAL
 #else  #else
     if (fx != fx){      if (fx != fx){
         xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */                          xxs=xxs/scale; /* Trying a smaller xx, closer to initial ax=0 */
         printf("|");                          printf("|");
         fprintf(ficlog,"|");                          fprintf(ficlog,"|");
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
         printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n",  axs, xxs, fx,fb, fa, xx, ax, bx);                          printf("\nLinmin NAN : input [axs=%lf:xxs=%lf], mnbrak outputs fx=%lf <(fb=%lf and fa=%lf) with xx=%lf in [ax=%lf:bx=%lf] \n",  axs, xxs, fx,fb, fa, xx, ax, bx);
 #endif  #endif
     }      }
   }while(fx != fx);    }while(fx != fx && xxs > 1.e-5);
 #endif  #endif
       
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
   printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);    printf("\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
   fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);    fprintf(ficlog,"\nLinmin after mnbrak: ax=%12.7f xx=%12.7f bx=%12.7f fa=%12.2f fx=%12.2f fb=%12.2f\n",  ax,xx,bx,fa,fx,fb);
 #endif  #endif
   #ifdef LINMINORIGINAL
   #else
           if(fb == fx){ /* Flat function in the direction */
                   xmin=xx;
       *flat=1;
           }else{
       *flat=0;
   #endif
   *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/    *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); /* Giving a bracketting triplet (ax, xx, bx), find a minimum, xmin, according to f1dim, *fret(xmin),*/
   /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */    /* fa = f(p[j] + ax * xi[j]), fx = f(p[j] + xx * xi[j]), fb = f(p[j] + bx * xi[j]) */
   /* fmin = f(p[j] + xmin * xi[j]) */    /* fmin = f(p[j] + xmin * xi[j]) */
   /* P+lambda n in that direction (lambdamin), with TOL between abscisses */    /* P+lambda n in that direction (lambdamin), with TOL between abscisses */
   /* f1dim(xmin): for (j=1;j<=ncom;j++) xt[j]=pcom[j]+xmin*xicom[j]; */    /* f1dim(xmin): for (j=1;j<=ncom;j++) xt[j]=pcom[j]+xmin*xicom[j]; */
 #ifdef DEBUG  #ifdef DEBUG
   printf("retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);    printf("retour brent from bracket (a=%lf fa=%lf, xx=%lf fx=%lf, b=%lf fb=%lf): fret=%lf xmin=%lf\n",ax,fa,xx,fx,bx,fb,*fret,xmin);
   fprintf(ficlog,"retour brent fret=%.12e xmin=%.12e\n",*fret,xmin);    fprintf(ficlog,"retour brent from bracket (a=%lf fa=%lf, xx=%lf fx=%lf, b=%lf fb=%lf): fret=%lf xmin=%lf\n",ax,fa,xx,fx,bx,fb,*fret,xmin);
   #endif
   #ifdef LINMINORIGINAL
   #else
                           }
 #endif  #endif
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
   printf("linmin end ");    printf("linmin end ");
Line 1835  such that failure to decrease by more th Line 1866  such that failure to decrease by more th
 output, p is set to the best point found, xi is the then-current direction set, fret is the returned  output, p is set to the best point found, xi is the then-current direction set, fret is the returned
 function value at p , and iter is the number of iterations taken. The routine linmin is used.  function value at p , and iter is the number of iterations taken. The routine linmin is used.
  */   */
   #ifdef LINMINORIGINAL
   #else
           int *flatdir; /* Function is vanishing in that direction */
           int flat=0; /* Function is vanishing in that direction */
   #endif
 void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret,   void powell(double p[], double **xi, int n, double ftol, int *iter, double *fret, 
             double (*func)(double []))               double (*func)(double [])) 
 {   { 
   void linmin(double p[], double xi[], int n, double *fret,   #ifdef LINMINORIGINAL
    void linmin(double p[], double xi[], int n, double *fret, 
               double (*func)(double []));                 double (*func)(double [])); 
   #else 
    void linmin(double p[], double xi[], int n, double *fret, 
                                                    double (*func)(double []),int *flat); 
   #endif
   int i,ibig,j;     int i,ibig,j; 
   double del,t,*pt,*ptt,*xit;    double del,t,*pt,*ptt,*xit;
   double directest;    double directest;
   double fp,fptt;    double fp,fptt;
   double *xits;    double *xits;
   int niterf, itmp;    int niterf, itmp;
   #ifdef LINMINORIGINAL
   #else
   
     flatdir=ivector(1,n); 
     for (j=1;j<=n;j++) flatdir[j]=0; 
   #endif
   
   pt=vector(1,n);     pt=vector(1,n); 
   ptt=vector(1,n);     ptt=vector(1,n); 
Line 1879  void powell(double p[], double **xi, int Line 1926  void powell(double p[], double **xi, int
       rforecast_time=rcurr_time;         rforecast_time=rcurr_time; 
       itmp = strlen(strcurr);        itmp = strlen(strcurr);
       if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */        if(strcurr[itmp-1]=='\n')  /* Windows outputs with a new line */
         strcurr[itmp-1]='\0';                                  strcurr[itmp-1]='\0';
       printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);        printf("\nConsidering the time needed for the last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);        fprintf(ficlog,"\nConsidering the time needed for this last iteration #%d: %ld seconds,\n",*iter,rcurr_time-rlast_time);
       for(niterf=10;niterf<=30;niterf+=10){        for(niterf=10;niterf<=30;niterf+=10){
         rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);                                  rforecast_time=rcurr_time+(niterf-*iter)*(rcurr_time-rlast_time);
         forecast_time = *localtime(&rforecast_time);                                  forecast_time = *localtime(&rforecast_time);
         strcpy(strfor,asctime(&forecast_time));                                  strcpy(strfor,asctime(&forecast_time));
         itmp = strlen(strfor);                                  itmp = strlen(strfor);
         if(strfor[itmp-1]=='\n')                                  if(strfor[itmp-1]=='\n')
         strfor[itmp-1]='\0';                                          strfor[itmp-1]='\0';
         printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);                                  printf("   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
         fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);                                  fprintf(ficlog,"   - if your program needs %d iterations to converge, convergence will be \n   reached in %s i.e.\n   on %s (current time is %s);\n",niterf, asc_diff_time(rforecast_time-rcurr_time,tmpout),strfor,strcurr);
       }        }
     }      }
     for (i=1;i<=n;i++) { /* For each direction i */      for (i=1;i<=n;i++) { /* For each direction i */
Line 1902  void powell(double p[], double **xi, int Line 1949  void powell(double p[], double **xi, int
 #endif  #endif
       printf("%d",i);fflush(stdout); /* print direction (parameter) i */        printf("%d",i);fflush(stdout); /* print direction (parameter) i */
       fprintf(ficlog,"%d",i);fflush(ficlog);        fprintf(ficlog,"%d",i);fflush(ficlog);
   #ifdef LINMINORIGINAL
       linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/        linmin(p,xit,n,fret,func); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
                                     /* Outputs are fret(new point p) p is updated and xit rescaled */  #else
         linmin(p,xit,n,fret,func,&flat); /* Point p[n]. xit[n] has been loaded for direction i as input.*/
                           flatdir[i]=flat; /* Function is vanishing in that direction i */
   #endif
                           /* Outputs are fret(new point p) p is updated and xit rescaled */
       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. */
         /* Unless the n directions are conjugate some gain in the determinant may be obtained */                                  /* Unless the n directions are conjugate some gain in the determinant may be obtained */
         /* with the new direction. */                                  /* with the new direction. */
         del=fabs(fptt-(*fret));                                   del=fabs(fptt-(*fret)); 
         ibig=i;                                   ibig=i; 
       }         } 
 #ifdef DEBUG  #ifdef DEBUG
       printf("%d %.12e",i,(*fret));        printf("%d %.12e",i,(*fret));
       fprintf(ficlog,"%d %.12e",i,(*fret));        fprintf(ficlog,"%d %.12e",i,(*fret));
       for (j=1;j<=n;j++) {        for (j=1;j<=n;j++) {
         xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);                                  xits[j]=FMAX(fabs(p[j]-pt[j]),1.e-5);
         printf(" x(%d)=%.12e",j,xit[j]);                                  printf(" x(%d)=%.12e",j,xit[j]);
         fprintf(ficlog," x(%d)=%.12e",j,xit[j]);                                  fprintf(ficlog," x(%d)=%.12e",j,xit[j]);
       }        }
       for(j=1;j<=n;j++) {        for(j=1;j<=n;j++) {
         printf(" p(%d)=%.12e",j,p[j]);                          printf(" p(%d)=%lf ",j,p[j]);
         fprintf(ficlog," p(%d)=%.12e",j,p[j]);                          fprintf(ficlog," p(%d)=%lf ",j,p[j]);
       }        }
       printf("\n");        printf("\n");
       fprintf(ficlog,"\n");        fprintf(ficlog,"\n");
Line 1931  void powell(double p[], double **xi, int Line 1983  void powell(double p[], double **xi, int
     /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */       /* Convergence test will use last linmin estimation (fret) and compare former iteration (fp) */ 
     /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */      /* But p and xit have been updated at the end of linmin, *fret corresponds to new p, xit  */
     /* New value of last point Pn is not computed, P(n-1) */      /* New value of last point Pn is not computed, P(n-1) */
         for(j=1;j<=n;j++) {
                             printf(" p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
                             fprintf(ficlog," p(%d)=%lf flat=%d ",j,p[j],flatdir[j]);
         }
         printf("\n");
         fprintf(ficlog,"\n");
   
     if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */      if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) { /* Did we reach enough precision? */
       /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */        /* We could compare with a chi^2. chisquare(0.95,ddl=1)=3.84 */
       /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */        /* By adding age*age in a model, the new -2LL should be lower and the difference follows a */
Line 1939  void powell(double p[], double **xi, int Line 1998  void powell(double p[], double **xi, int
       /* By adding age*age and V1*age the gain (-2LL) should be more than 5.99 (ddl=2) */        /* By adding age*age and V1*age the gain (-2LL) should be more than 5.99 (ddl=2) */
       /* By using V1+V2+V3, the gain should be  7.82, compared with basic 1+age. */        /* By using V1+V2+V3, the gain should be  7.82, compared with basic 1+age. */
       /* By adding 10 parameters more the gain should be 18.31 */        /* By adding 10 parameters more the gain should be 18.31 */
                           
       /* Starting the program with initial values given by a former maximization will simply change */        /* Starting the program with initial values given by a former maximization will simply change */
       /* the scales of the directions and the directions, because the are reset to canonical directions */        /* the scales of the directions and the directions, because the are reset to canonical directions */
       /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */        /* Thus the first calls to linmin will give new points and better maximizations until fp-(*fret) is */
Line 1967  void powell(double p[], double **xi, int Line 2026  void powell(double p[], double **xi, int
       }        }
 #endif  #endif
   
   #ifdef LINMINORIGINAL
   #else
         free_ivector(flatdir,1,n); 
   #endif
       free_vector(xit,1,n);         free_vector(xit,1,n); 
       free_vector(xits,1,n);         free_vector(xits,1,n); 
       free_vector(ptt,1,n);         free_vector(ptt,1,n); 
Line 1981  void powell(double p[], double **xi, int Line 2043  void powell(double p[], double **xi, int
       pt[j]=p[j];         pt[j]=p[j]; 
     }       } 
     fptt=(*func)(ptt); /* f_3 */      fptt=(*func)(ptt); /* f_3 */
 #ifdef POWELLF1F3  #ifdef POWELLNOF3INFF1TEST    /* skips test F3 <F1 */
 #else  #else
     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 */
 #endif  #endif
Line 1990  void powell(double p[], double **xi, int Line 2052  void powell(double p[], double **xi, int
       /* 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 */
       /* 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 or directest <0 */
         /* also  lamda^2=(f1-f2)^2/mu² is a parasite solution of powell */
         /* For powell, inclusion of this average direction is only if t(del)<0 or del inbetween mu^2 and lambda^2 */
       /* 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); */
         /*  Even if f3 <f1, directest can be negative and t >0 */
         /* mu² and del² are equal when f3=f1 */
                           /* f3 < f1 : mu² < del <= lambda^2 both test are equivalent */
                           /* f3 < f1 : mu² < lambda^2 < del then directtest is negative and powell t is positive */
                           /* f3 > f1 : lambda² < mu^2 < del then t is negative and directest >0  */
                           /* f3 > f1 : lambda² < del < mu^2 then t is positive and directest >0  */
 #ifdef NRCORIGINAL  #ifdef NRCORIGINAL
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)- del*SQR(fp-fptt); /* Original Numerical Recipes in C*/
 #else  #else
Line 2013  void powell(double p[], double **xi, int Line 2083  void powell(double p[], double **xi, int
       if (t < 0.0) { /* Then we use it for new direction */        if (t < 0.0) { /* Then we use it for new direction */
 #else  #else
       if (directest*t < 0.0) { /* Contradiction between both tests */        if (directest*t < 0.0) { /* Contradiction between both tests */
         printf("directest= %.12lf (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt,del);                                  printf("directest= %.12lf (if <0 we include P0 Pn as new direction), 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 (if <0 we include P0 Pn as new direction), t= %.12lf, f1= %.12lf,f2= %.12lf,f3= %.12lf, del= %.12lf\n",directest, t, fp,(*fret),fptt, del);          fprintf(ficlog,"directest= %.12lf (if directest<0 or t<0 we include P0 Pn as new direction), 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 */
 #endif  #endif
 #ifdef DEBUGLINMIN  #ifdef DEBUGLINMIN
         printf("Before linmin in direction P%d-P0\n",n);                                  printf("Before linmin in direction P%d-P0\n",n);
         for (j=1;j<=n;j++) {                                   for (j=1;j<=n;j++) {
           printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);                                          printf(" Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
           fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);                                          fprintf(ficlog," Before xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
           if(j % ncovmodel == 0){                                          if(j % ncovmodel == 0){
             printf("\n");                                                  printf("\n");
             fprintf(ficlog,"\n");                                                  fprintf(ficlog,"\n");
           }                                          }
         }                                  }
 #endif  #endif
         linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/  #ifdef LINMINORIGINAL
 #ifdef DEBUGLINMIN                                  linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
         for (j=1;j<=n;j++) {   #else
           printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);                                  linmin(p,xit,n,fret,func,&flat); /* computes minimum on the extrapolated direction: changes p and rescales xit.*/
           fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);                                  flatdir[i]=flat; /* Function is vanishing in that direction i */
           if(j % ncovmodel == 0){  
             printf("\n");  
             fprintf(ficlog,"\n");  
           }  
         }  
 #endif  #endif
         for (j=1;j<=n;j++) {   
           xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */  
           xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */  
         }  
         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 new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);  
   
   #ifdef DEBUGLINMIN
                                   for (j=1;j<=n;j++) { 
                                           printf("After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
                                           fprintf(ficlog,"After xit[%d]= %12.7f p[%d]= %12.7f",j,xit[j],j,p[j]);
                                           if(j % ncovmodel == 0){
                                                   printf("\n");
                                                   fprintf(ficlog,"\n");
                                           }
                                   }
   #endif
                                   for (j=1;j<=n;j++) { 
                                           xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
                                           xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
                                   }
   #ifdef LINMINORIGINAL
   #else
                                   printf("Flat directions\n");
                                   fprintf(ficlog,"Flat directions\n");
                                   for (j=1;j<=n;j++) { 
                                           printf("flatdir[%d]=%d ",j,flatdir[j]);
                                           fprintf(ficlog,"flatdir[%d]=%d ",j,flatdir[j]);
           }
                                   printf("\n");
                                   fprintf(ficlog,"\n");
   #endif
                                   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 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);
         fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);                                  fprintf(ficlog,"Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
         for(j=1;j<=n;j++){                                  for(j=1;j<=n;j++){
           printf(" %.12e",xit[j]);                                          printf(" %lf",xit[j]);
           fprintf(ficlog," %.12e",xit[j]);                                          fprintf(ficlog," %lf",xit[j]);
         }                                  }
         printf("\n");                                  printf("\n");
         fprintf(ficlog,"\n");                                  fprintf(ficlog,"\n");
 #endif  #endif
       } /* end of t or directest negative */        } /* end of t or directest negative */
 #ifdef POWELLF1F3  #ifdef POWELLNOF3INFF1TEST
 #else  #else
     } /* end if (fptt < fp)  */      } /* end if (fptt < fp)  */
 #endif  #endif
Line 2337  double **pmij(double **ps, double *cov, Line 2424  double **pmij(double **ps, double *cov,
   /*double t34;*/    /*double t34;*/
   int i,j, nc, ii, jj;    int i,j, nc, ii, jj;
   
         for(i=1; i<= nlstate; i++){    for(i=1; i<= nlstate; i++){
                 for(j=1; j<i;j++){      for(j=1; j<i;j++){
                         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){        for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
                                 /*lnpijopii += param[i][j][nc]*cov[nc];*/          /*lnpijopii += param[i][j][nc]*cov[nc];*/
                                 lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];          lnpijopii += x[nc+((i-1)*(nlstate+ndeath-1)+j-1)*ncovmodel]*cov[nc];
                                 /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */          /*       printf("Int j<i s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
                         }        }
                         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */        ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
                         /*      printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */        /*        printf("s1=%.17e, lnpijopii=%.17e\n",s1,lnpijopii); */
                 }      }
                 for(j=i+1; j<=nlstate+ndeath;j++){      for(j=i+1; j<=nlstate+ndeath;j++){
                         for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){        for (nc=1, lnpijopii=0.;nc <=ncovmodel; nc++){
                                 /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/          /*lnpijopii += x[(i-1)*nlstate*ncovmodel+(j-2)*ncovmodel+nc+(i-1)*(ndeath-1)*ncovmodel]*cov[nc];*/
                                 lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];          lnpijopii += x[nc + ((i-1)*(nlstate+ndeath-1)+(j-2))*ncovmodel]*cov[nc];
                                 /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */          /*        printf("Int j>i s1=%.17e, lnpijopii=%.17e %lx %lx\n",s1,lnpijopii,s1,lnpijopii); */
                         }        }
                         ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */        ps[i][j]=lnpijopii; /* In fact ln(pij/pii) */
                 }      }
         }    }
       
         for(i=1; i<= nlstate; i++){    for(i=1; i<= nlstate; i++){
                 s1=0;      s1=0;
                 for(j=1; j<i; j++){      for(j=1; j<i; j++){
                         s1+=exp(ps[i][j]); /* In fact sums pij/pii */        s1+=exp(ps[i][j]); /* In fact sums pij/pii */
                         /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */        /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
                 }      }
                 for(j=i+1; j<=nlstate+ndeath; j++){      for(j=i+1; j<=nlstate+ndeath; j++){
                         s1+=exp(ps[i][j]); /* In fact sums pij/pii */        s1+=exp(ps[i][j]); /* In fact sums pij/pii */
                         /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */        /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
                 }      }
                 /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */      /* s1= sum_{j<>i} pij/pii=(1-pii)/pii and thus pii is known from s1 */
                 ps[i][i]=1./(s1+1.);      ps[i][i]=1./(s1+1.);
                 /* Computing other pijs */      /* Computing other pijs */
                 for(j=1; j<i; j++)      for(j=1; j<i; j++)
                         ps[i][j]= exp(ps[i][j])*ps[i][i];        ps[i][j]= exp(ps[i][j])*ps[i][i];
                 for(j=i+1; j<=nlstate+ndeath; j++)      for(j=i+1; j<=nlstate+ndeath; j++)
                         ps[i][j]= exp(ps[i][j])*ps[i][i];        ps[i][j]= exp(ps[i][j])*ps[i][i];
                 /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */      /* ps[i][nlstate+1]=1.-s1- ps[i][i];*/ /* Sum should be 1 */
         } /* end i */    } /* end i */
       
         for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){    for(ii=nlstate+1; ii<= nlstate+ndeath; ii++){
                 for(jj=1; jj<= nlstate+ndeath; jj++){      for(jj=1; jj<= nlstate+ndeath; jj++){
                         ps[ii][jj]=0;        ps[ii][jj]=0;
                         ps[ii][ii]=1;        ps[ii][ii]=1;
                 }      }
         }    }
       
       
         /* for(ii=1; ii<= nlstate+ndeath; ii++){ */    /* for(ii=1; ii<= nlstate+ndeath; ii++){ */
         /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */    /*   for(jj=1; jj<= nlstate+ndeath; jj++){ */
         /*      printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */    /*    printf(" pmij  ps[%d][%d]=%lf ",ii,jj,ps[ii][jj]); */
         /*   } */    /*   } */
         /*   printf("\n "); */    /*   printf("\n "); */
         /* } */    /* } */
         /* printf("\n ");printf("%lf ",cov[2]);*/    /* printf("\n ");printf("%lf ",cov[2]);*/
         /*    /*
                 for(i=1; i<= npar; i++) printf("%f ",x[i]);      for(i=1; i<= npar; i++) printf("%f ",x[i]);
                 goto end;*/                  goto end;*/
         return ps;    return ps;
 }  }
   
 /*************** backward transition probabilities ***************/   /*************** backward transition probabilities ***************/ 
Line 2790  double func( double *x) Line 2877  double func( double *x)
     printf(" %d\n",s[4][i]);      printf(" %d\n",s[4][i]);
   */    */
   
   ++countcallfunc;          ++countcallfunc;
   
   cov[1]=1.;          cov[1]=1.;
   
   for(k=1; k<=nlstate; k++) ll[k]=0.;    for(k=1; k<=nlstate; k++) ll[k]=0.;
   
Line 2894  double func( double *x) Line 2981  double func( double *x)
 /*        else */  /*        else */
 /*          lli=log(out[s1][s2] - savm[s1][s2]); */  /*          lli=log(out[s1][s2] - savm[s1][s2]); */
 /* #endif */  /* #endif */
           lli=log(out[s1][s2] - savm[s1][s2]);                                          lli=log(out[s1][s2] - savm[s1][s2]);
                       
         } else if  ( s2==-1 ) { /* alive */                                  } else if  ( s2==-1 ) { /* alive */
           for (j=1,survp=0. ; j<=nlstate; j++)                                           for (j=1,survp=0. ; j<=nlstate; j++) 
             survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];                                                  survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
           /*survp += out[s1][j]; */                                          /*survp += out[s1][j]; */
           lli= log(survp);                                          lli= log(survp);
         }                                  }
         else if  (s2==-4) {                                   else if  (s2==-4) { 
           for (j=3,survp=0. ; j<=nlstate; j++)                                            for (j=3,survp=0. ; j<=nlstate; j++)  
             survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];                                                  survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
           lli= log(survp);                                           lli= log(survp); 
         }                                   } 
         else if  (s2==-5) {                                   else if  (s2==-5) { 
           for (j=1,survp=0. ; j<=2; j++)                                            for (j=1,survp=0. ; j<=2; j++)  
             survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];                                                  survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
           lli= log(survp);                                           lli= log(survp); 
         }                                   } 
         else{                                  else{
           lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */                                          lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
           /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */                                          /*  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2]));*/ /* linear interpolation */
         }                                   } 
         /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/                                  /*lli=(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]);*/
         /*if(lli ==000.0)*/                                  /*if(lli ==000.0)*/
         /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */                                  /*printf("bbh= %f lli=%f savm=%f out=%f %d\n",bbh,lli,savm[s1][s2], out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]],i); */
         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)){ */                                  /* 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); */                                  /*   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]); */                                  /*   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){
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){                  for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];                          for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){                          for(mi=1; mi<= wav[i]-1; mi++){
         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);
             savm[ii][j]=(ii==j ? 1.0 : 0.0);                                                  savm[ii][j]=(ii==j ? 1.0 : 0.0);
           }                                          }
         for(d=0; d<=dh[mi][i]; d++){                                  for(d=0; d<=dh[mi][i]; d++){
           newm=savm;                                          newm=savm;
           agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;                                          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
           cov[2]=agexact;                                          cov[2]=agexact;
           if(nagesqr==1)                                          if(nagesqr==1)
             cov[3]= agexact*agexact;                                                  cov[3]= agexact*agexact;
           for (kk=1; kk<=cptcovage;kk++) {                                          for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;                                                  cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
           }                                          }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,                                          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                                                                                           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           savm=oldm;                                          savm=oldm;
           oldm=newm;                                          oldm=newm;
         } /* end mult */                                  } /* end mult */
               
         s1=s[mw[mi][i]][i];                                  s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];                                  s2=s[mw[mi+1][i]][i];
         bbh=(double)bh[mi][i]/(double)stepm;                                   bbh=(double)bh[mi][i]/(double)stepm; 
         lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */                                  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*(savm[s1][s2])):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
         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;
       } /* end of wave */                          } /* end of wave */
     } /* end of individual */                  } /* end of individual */
   }  else if(mle==3){  /* exponential inter-extrapolation */          }  else if(mle==3){  /* exponential inter-extrapolation */
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){                  for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];                          for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){                          for(mi=1; mi<= wav[i]-1; mi++){
         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);
             savm[ii][j]=(ii==j ? 1.0 : 0.0);                                                  savm[ii][j]=(ii==j ? 1.0 : 0.0);
           }                                          }
         for(d=0; d<dh[mi][i]; d++){                                  for(d=0; d<dh[mi][i]; d++){
           newm=savm;                                          newm=savm;
           agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;                                          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
           cov[2]=agexact;                                          cov[2]=agexact;
           if(nagesqr==1)                                          if(nagesqr==1)
             cov[3]= agexact*agexact;                                                  cov[3]= agexact*agexact;
           for (kk=1; kk<=cptcovage;kk++) {                                          for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;                                                  cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
           }                                          }
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,                                          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                                                                                           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           savm=oldm;                                          savm=oldm;
           oldm=newm;                                          oldm=newm;
         } /* end mult */                                  } /* end mult */
               
         s1=s[mw[mi][i]][i];                                  s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];                                  s2=s[mw[mi+1][i]][i];
         bbh=(double)bh[mi][i]/(double)stepm;                                   bbh=(double)bh[mi][i]/(double)stepm; 
         lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */                                  lli= (savm[s1][s2]>1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
         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;
       } /* end of wave */                          } /* end of wave */
     } /* end of individual */                  } /* end of individual */
   }else if (mle==4){  /* ml=4 no inter-extrapolation */          }else if (mle==4){  /* ml=4 no inter-extrapolation */
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){                  for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];                          for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){                          for(mi=1; mi<= wav[i]-1; mi++){
         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);
             savm[ii][j]=(ii==j ? 1.0 : 0.0);                                                  savm[ii][j]=(ii==j ? 1.0 : 0.0);
           }                                          }
         for(d=0; d<dh[mi][i]; d++){                                  for(d=0; d<dh[mi][i]; d++){
           newm=savm;                                          newm=savm;
           agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;                                          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
           cov[2]=agexact;                                          cov[2]=agexact;
           if(nagesqr==1)                                          if(nagesqr==1)
             cov[3]= agexact*agexact;                                                  cov[3]= agexact*agexact;
           for (kk=1; kk<=cptcovage;kk++) {                                          for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;                                                  cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
           }                                          }
                   
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,                                          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                                                                                           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           savm=oldm;                                          savm=oldm;
           oldm=newm;                                          oldm=newm;
         } /* end mult */                                  } /* end mult */
               
         s1=s[mw[mi][i]][i];                                  s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];                                  s2=s[mw[mi+1][i]][i];
         if( s2 > nlstate){                                   if( s2 > nlstate){ 
           lli=log(out[s1][s2] - savm[s1][s2]);                                          lli=log(out[s1][s2] - savm[s1][s2]);
         } else if  ( s2==-1 ) { /* alive */                                  } else if  ( s2==-1 ) { /* alive */
           for (j=1,survp=0. ; j<=nlstate; j++)                                           for (j=1,survp=0. ; j<=nlstate; j++) 
             survp += out[s1][j];                                                  survp += out[s1][j];
           lli= log(survp);                                          lli= log(survp);
         }else{                                  }else{
           lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */                                          lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
         }                                  }
         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;
 /*      printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */  /*      printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       } /* end of wave */                          } /* end of wave */
     } /* end of individual */                  } /* end of individual */
   }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */          }else{  /* ml=5 no inter-extrapolation no jackson =0.8a */
     for (i=1,ipmx=0, sw=0.; i<=imx; i++){                  for (i=1,ipmx=0, sw=0.; i<=imx; i++){
       for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];                          for (k=1; k<=cptcovn;k++) cov[2+nagesqr+k]=covar[Tvar[k]][i];
       for(mi=1; mi<= wav[i]-1; mi++){                          for(mi=1; mi<= wav[i]-1; mi++){
         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);
             savm[ii][j]=(ii==j ? 1.0 : 0.0);                                                  savm[ii][j]=(ii==j ? 1.0 : 0.0);
           }                                          }
         for(d=0; d<dh[mi][i]; d++){                                  for(d=0; d<dh[mi][i]; d++){
           newm=savm;                                          newm=savm;
           agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;                                          agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
           cov[2]=agexact;                                          cov[2]=agexact;
           if(nagesqr==1)                                          if(nagesqr==1)
             cov[3]= agexact*agexact;                                                  cov[3]= agexact*agexact;
           for (kk=1; kk<=cptcovage;kk++) {                                          for (kk=1; kk<=cptcovage;kk++) {
             cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;                                                  cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
           }                                          }
                   
           out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,                                          out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                        1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                                                                                           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
           savm=oldm;                                          savm=oldm;
           oldm=newm;                                          oldm=newm;
         } /* end mult */                                  } /* end mult */
               
         s1=s[mw[mi][i]][i];                                  s1=s[mw[mi][i]][i];
         s2=s[mw[mi+1][i]][i];                                  s2=s[mw[mi+1][i]][i];
         lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */                                  lli=log(out[s[mw[mi][i]][i]][s[mw[mi+1][i]][i]]); /* Original formula */
         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;
         /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/                                  /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);*/
       } /* end of wave */                          } /* end of wave */
     } /* end of individual */                  } /* end of individual */
   } /* End of if */          } /* End of if */
   for(k=1,l=0.; k<=nlstate; k++) l += ll[k];          for(k=1,l=0.; k<=nlstate; k++) l += ll[k];
   /* printf("l1=%f l2=%f ",ll[1],ll[2]); */          /* printf("l1=%f l2=%f ",ll[1],ll[2]); */
   l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */          l= l*ipmx/sw; /* To get the same order of magnitude as if weight=1 for every body */
   return -l;          return -l;
 }  }
   
 /*************** log-likelihood *************/  /*************** log-likelihood *************/
Line 3112  double funcone( double *x) Line 3199  double funcone( double *x)
       agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */        agebegin=agev[mw[mi][i]][i]; /* Age at beginning of effective wave */
       ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */        ageend=agev[mw[mi][i]][i] + (dh[mi][i])*stepm/YEARM; /* Age at end of effective wave and at the end of transition */
       for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */        for(d=0; d<dh[mi][i]; d++){  /* Delay between two effective waves */
         /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]                                  /*dh[m][i] or dh[mw[mi][i]][i] is the delay between two effective waves m=mw[mi][i]
           and mw[mi+1][i]. dh depends on stepm.*/                                          and mw[mi+1][i]. dh depends on stepm.*/
         newm=savm;                                  newm=savm;
         agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;                                  agexact=agev[mw[mi][i]][i]+d*stepm/YEARM;
         cov[2]=agexact;                                  cov[2]=agexact;
         if(nagesqr==1)                                  if(nagesqr==1)
           cov[3]= agexact*agexact;                                          cov[3]= agexact*agexact;
         for (kk=1; kk<=cptcovage;kk++) {                                  for (kk=1; kk<=cptcovage;kk++) {
           cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;                                          cov[Tage[kk]+2+nagesqr]=covar[Tvar[Tage[kk]]][i]*agexact;
         }                                  }
                                   
         /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */                                  /* savm=pmij(pmmij,cov,ncovmodel,x,nlstate); */
         out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,                                  out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,
                      1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));                                                                                   1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate));
         /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */                                  /* out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath, */
         /*           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */                                  /*           1,nlstate+ndeath,pmij(pmmij,cov,ncovmodel,x,nlstate)); */
         savm=oldm;                                  savm=oldm;
         oldm=newm;                                  oldm=newm;
       } /* end mult */        } /* end mult */
               
       s1=s[mw[mi][i]][i];        s1=s[mw[mi][i]][i];
Line 3143  double funcone( double *x) Line 3230  double funcone( double *x)
        * is higher than the multiple of stepm and negative otherwise.         * is higher than the multiple of stepm and negative otherwise.
        */         */
       if( s2 > nlstate && (mle <5) ){  /* Jackson */        if( s2 > nlstate && (mle <5) ){  /* Jackson */
         lli=log(out[s1][s2] - savm[s1][s2]);                                  lli=log(out[s1][s2] - savm[s1][s2]);
       } else if  ( s2==-1 ) { /* alive */        } else if  ( s2==-1 ) { /* alive */
         for (j=1,survp=0. ; j<=nlstate; j++)                                   for (j=1,survp=0. ; j<=nlstate; j++) 
           survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];                                          survp += (1.+bbh)*out[s1][j]- bbh*savm[s1][j];
         lli= log(survp);                                  lli= log(survp);
       }else if (mle==1){        }else if (mle==1){
         lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */                                  lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
       } else if(mle==2){        } else if(mle==2){
         lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */                                  lli= (savm[s1][s2]>(double)1.e-8 ?log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* linear interpolation */
       } else if(mle==3){  /* exponential inter-extrapolation */        } else if(mle==3){  /* exponential inter-extrapolation */
         lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */                                  lli= (savm[s1][s2]>(double)1.e-8 ?(1.+bbh)*log(out[s1][s2])- bbh*log(savm[s1][s2]):log((1.+bbh)*out[s1][s2])); /* exponential inter-extrapolation */
       } else if (mle==4){  /* mle=4 no inter-extrapolation */        } else if (mle==4){  /* mle=4 no inter-extrapolation */
         lli=log(out[s1][s2]); /* Original formula */                                  lli=log(out[s1][s2]); /* Original formula */
       } else{  /* mle=0 back to 1 */        } else{  /* mle=0 back to 1 */
         lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */                                  lli= log((1.+bbh)*out[s1][s2]- bbh*savm[s1][s2]); /* linear interpolation */
         /*lli=log(out[s1][s2]); */ /* Original formula */                                  /*lli=log(out[s1][s2]); */ /* Original formula */
       } /* End of if */        } /* End of if */
       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;
       /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */        /*printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
       if(globpr){        if(globpr){
         fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\                                  fprintf(ficresilk,"%9ld %6.1f %6.1f %6d %2d %2d %2d %2d %3d %11.6f %8.4f %8.3f\
  %11.6f %11.6f %11.6f ", \   %11.6f %11.6f %11.6f ", \
                 num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,                                                                  num[i], agebegin, ageend, i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],weight[i]*gipmx/gsw,
                 2*weight[i]*lli,out[s1][s2],savm[s1][s2]);                                                                  2*weight[i]*lli,out[s1][s2],savm[s1][s2]);
         for(k=1,llt=0.,l=0.; k<=nlstate; k++){                                  for(k=1,llt=0.,l=0.; k<=nlstate; k++){
           llt +=ll[k]*gipmx/gsw;                                          llt +=ll[k]*gipmx/gsw;
           fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);                                          fprintf(ficresilk," %10.6f",-ll[k]*gipmx/gsw);
         }                                  }
         fprintf(ficresilk," %10.6f\n", -llt);                                  fprintf(ficresilk," %10.6f\n", -llt);
       }        }
     } /* end of wave */      } /* end of wave */
   } /* end of individual */    } /* end of individual */
Line 3685  void pstamp(FILE *fichier) Line 3772  void pstamp(FILE *fichier)
                                                                          int firstpass,  int lastpass, int stepm, int weightopt, char model[])                                                                           int firstpass,  int lastpass, int stepm, int weightopt, char model[])
  {  /* Some frequencies */   {  /* Some frequencies */
       
          int i, m, jk, j1, bool, z1,j;     int i, m, jk, j1, bool, z1,j;
          int iind=0, iage=0;     int iind=0, iage=0;
          int mi; /* Effective wave */     int mi; /* Effective wave */
          int first;     int first;
          double ***freq; /* Frequencies */     double ***freq; /* Frequencies */
          double *pp, **prop, *posprop, *pospropt;     double *pp, **prop, *posprop, *pospropt;
          double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0;     double pos=0., posproptt=0., pospropta=0., k2, dateintsum=0,k2cpt=0;
          char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];     char fileresp[FILENAMELENGTH], fileresphtm[FILENAMELENGTH], fileresphtmfr[FILENAMELENGTH];
          double agebegin, ageend;     double agebegin, ageend;
           
          pp=vector(1,nlstate);     pp=vector(1,nlstate);
          prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE);      prop=matrix(1,nlstate,iagemin-AGEMARGE,iagemax+3+AGEMARGE); 
          posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */      posprop=vector(1,nlstate); /* Counting the number of transition starting from a live state per age */ 
          pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */      pospropt=vector(1,nlstate); /* Counting the number of transition starting from a live state */ 
          /* prop=matrix(1,nlstate,iagemin,iagemax+3); */     /* prop=matrix(1,nlstate,iagemin,iagemax+3); */
          strcpy(fileresp,"P_");     strcpy(fileresp,"P_");
          strcat(fileresp,fileresu);     strcat(fileresp,fileresu);
          /*strcat(fileresphtm,fileresu);*/     /*strcat(fileresphtm,fileresu);*/
          if((ficresp=fopen(fileresp,"w"))==NULL) {     if((ficresp=fopen(fileresp,"w"))==NULL) {
                  printf("Problem with prevalence resultfile: %s\n", fileresp);       printf("Problem with prevalence resultfile: %s\n", fileresp);
                  fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);       fprintf(ficlog,"Problem with prevalence resultfile: %s\n", fileresp);
                  exit(0);       exit(0);
          }     }
   
          strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));     strcpy(fileresphtm,subdirfext(optionfilefiname,"PHTM_",".htm"));
          if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {     if((ficresphtm=fopen(fileresphtm,"w"))==NULL) {
                  printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));       printf("Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
                  fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));       fprintf(ficlog,"Problem with prevalence HTM resultfile '%s' with errno='%s'\n",fileresphtm,strerror(errno));
                  fflush(ficlog);       fflush(ficlog);
                  exit(70);        exit(70); 
          }     }
          else{     else{
                  fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \       fprintf(ficresphtm,"<html><head>\n<title>IMaCh PHTM_ %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
                                                  fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);               fileresphtm,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
          }     }
          fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);     fprintf(ficresphtm,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies and prevalence by age at begin of transition</h4>\n",fileresphtm, fileresphtm);
           
          strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));     strcpy(fileresphtmfr,subdirfext(optionfilefiname,"PHTMFR_",".htm"));
          if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {     if((ficresphtmfr=fopen(fileresphtmfr,"w"))==NULL) {
                  printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));       printf("Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
                  fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));       fprintf(ficlog,"Problem with frequency table HTM resultfile '%s' with errno='%s'\n",fileresphtmfr,strerror(errno));
                  fflush(ficlog);       fflush(ficlog);
                  exit(70);        exit(70); 
          }     }
          else{     else{
                  fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \       fprintf(ficresphtmfr,"<html><head>\n<title>IMaCh PHTM_Frequency table %s</title></head>\n <body><font size=\"2\">%s <br> %s</font> \
 <hr size=\"2\" color=\"#EC5E5E\"> \n\  <hr size=\"2\" color=\"#EC5E5E\"> \n\
 Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\  Title=%s <br>Datafile=%s Firstpass=%d Lastpass=%d Stepm=%d Weight=%d Model=1+age+%s<br>\n",\
                                                  fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);               fileresphtmfr,version,fullversion,title,datafile,firstpass,lastpass,stepm, weightopt, model);
          }     }
          fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);     fprintf(ficresphtmfr,"Current page is file <a href=\"%s\">%s</a><br>\n\n<h4>Frequencies of all effective transitions by age at begin of transition </h4>Unknown status is -1<br/>\n",fileresphtmfr, fileresphtmfr);
   
          freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);     freq= ma3x(-5,nlstate+ndeath,-5,nlstate+ndeath,iagemin-AGEMARGE,iagemax+3+AGEMARGE);
          j1=0;     j1=0;
       
          j=cptcoveff;     j=cptcoveff;
          if (cptcovn<1) {j=1;ncodemax[1]=1;}     if (cptcovn<1) {j=1;ncodemax[1]=1;}
   
          first=1;     first=1;
   
          /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:     /* Detects if a combination j1 is empty: for a multinomial variable like 3 education levels:
                         reference=low_education V1=0,V2=0        reference=low_education V1=0,V2=0
                         med_educ                V1=1 V2=0,         med_educ                V1=1 V2=0, 
                         high_educ               V1=0 V2=1        high_educ               V1=0 V2=1
                         Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff         Then V1=1 and V2=1 is a noisy combination that we want to exclude for the list 2**cptcoveff 
          */     */
   
          for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */     for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){ /* Loop on covariates combination */
                  posproptt=0.;       posproptt=0.;
                  /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
                          scanf("%d", i);*/         scanf("%d", i);*/
                  for (i=-5; i<=nlstate+ndeath; i++)         for (i=-5; i<=nlstate+ndeath; i++)  
                          for (jk=-5; jk<=nlstate+ndeath; jk++)           for (jk=-5; jk<=nlstate+ndeath; jk++)  
                                  for(m=iagemin; m <= iagemax+3; m++)           for(m=iagemin; m <= iagemax+3; m++)
                                          freq[i][jk][m]=0;             freq[i][jk][m]=0;
               
                  for (i=1; i<=nlstate; i++)  {       for (i=1; i<=nlstate; i++)  {
                          for(m=iagemin; m <= iagemax+3; m++)         for(m=iagemin; m <= iagemax+3; m++)
                                  prop[i][m]=0;           prop[i][m]=0;
                          posprop[i]=0;         posprop[i]=0;
                          pospropt[i]=0;         pospropt[i]=0;
                  }       }
               
                  dateintsum=0;       dateintsum=0;
                  k2cpt=0;       k2cpt=0;
        /* For that comination of covariate j1, we count and print the frequencies */
                  for (iind=1; iind<=imx; iind++) { /* For each individual iind */       for (iind=1; iind<=imx; iind++) { /* For each individual iind */
                          bool=1;         bool=1;
                          if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */         if (cptcovn>0) { /* Filter is here: Must be looked at for model=V1+V2+V3+V4 */
                                  for (z1=1; z1<=cptcoveff; z1++) {                 for (z1=1; z1<=cptcoveff; z1++) {      
                                          if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){             if (covar[Tvaraff[z1]][iind]!= nbcode[Tvaraff[z1]][codtabm(j1,z1)]){
                                                  /* Tests if the value of each of the covariates of i is equal to filter j1 */               /* Tests if the value of each of the covariates of i is equal to filter j1 */
                                                  bool=0;               bool=0;
                                                  /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n",                /* printf("bool=%d i=%d, z1=%d, Tvaraff[%d]=%d, covar[Tvarff][%d]=%2f, codtabm(%d,%d)=%d, nbcode[Tvaraff][codtabm(%d,%d)=%d, j1=%d\n", 
                 bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),                  bool,i,z1, z1, Tvaraff[z1],i,covar[Tvaraff[z1]][i],j1,z1,codtabm(j1,z1),
                 j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/                  j1,z1,nbcode[Tvaraff[z1]][codtabm(j1,z1)],j1);*/
                                                  /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/               /* For j1=7 in V1+V2+V3+V4 = 0 1 1 0 and codtabm(7,3)=1 and nbcde[3][?]=1*/
                                          }              } 
                                  } /* end z1 */           } /* end z1 */
                          } /* cptcovn > 0 */         } /* cptcovn > 0 */
   
                          if (bool==1){         if (bool==1){ /* We selected an individual iin satisfying combination j1 */
                                  /* for(m=firstpass; m<=lastpass; m++){ */           /* for(m=firstpass; m<=lastpass; m++){ */
                                  for(mi=1; mi<wav[iind];mi++){           for(mi=1; mi<wav[iind];mi++){
                                          m=mw[mi][iind];             m=mw[mi][iind];
                                          /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]             /* dh[m][iind] or dh[mw[mi][iind]][iind] is the delay between two effective (mi) waves m=mw[mi][iind]
                                                         and mw[mi+1][iind]. dh depends on stepm. */                and mw[mi+1][iind]. dh depends on stepm. */
                                          agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/             agebegin=agev[m][iind]; /* Age at beginning of wave before transition*/
                                          ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */             ageend=agev[m][iind]+(dh[m][iind])*stepm/YEARM; /* Age at end of wave and transition */
                                          if(m >=firstpass && m <=lastpass){             if(m >=firstpass && m <=lastpass){
                                                  k2=anint[m][iind]+(mint[m][iind]/12.);               k2=anint[m][iind]+(mint[m][iind]/12.);
                                                  /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/               /*if ((k2>=dateprev1) && (k2<=dateprev2)) {*/
                                                  if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */               if(agev[m][iind]==0) agev[m][iind]=iagemax+1;  /* All ages equal to 0 are in iagemax+1 */
                                                  if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */               if(agev[m][iind]==1) agev[m][iind]=iagemax+2;  /* All ages equal to 1 are in iagemax+2 */
                                                  if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */               if (s[m][iind]>0 && s[m][iind]<=nlstate)  /* If status at wave m is known and a live state */
                                                          prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */                 prop[s[m][iind]][(int)agev[m][iind]] += weight[iind];  /* At age of beginning of transition, where status is known */
                                                  if (m<lastpass) {               if (m<lastpass) {
                                                          /* if(s[m][iind]==4 && s[m+1][iind]==4) */                 /* if(s[m][iind]==4 && s[m+1][iind]==4) */
                                                          /*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */                 /*   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind]); */
                                                          if(s[m][iind]==-1)                 if(s[m][iind]==-1)
                                                                  printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));                   printf(" num=%ld m=%d, iind=%d s1=%d s2=%d agev at m=%d agebegin=%.2f ageend=%.2f, agemed=%d\n", num[iind], m, iind,s[m][iind],s[m+1][iind], (int)agev[m][iind],agebegin, ageend, (int)((agebegin+ageend)/2.));
                                                          freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */                 freq[s[m][iind]][s[m+1][iind]][(int)agev[m][iind]] += weight[iind]; /* At age of beginning of transition, where status is known */
                                                          /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */                 /* freq[s[m][iind]][s[m+1][iind]][(int)((agebegin+ageend)/2.)] += weight[iind]; */
                                                          freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */                 freq[s[m][iind]][s[m+1][iind]][iagemax+3] += weight[iind]; /* Total is in iagemax+3 *//* At age of beginning of transition, where status is known */
                                                  }               }
                                          }               }  
                                          if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {             if ((agev[m][iind]>1) && (agev[m][iind]< (iagemax+3)) && (anint[m][iind]!=9999) && (mint[m][iind]!=99)) {
                                                  dateintsum=dateintsum+k2;               dateintsum=dateintsum+k2;
                                                  k2cpt++;               k2cpt++;
                                                  /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */               /* printf("iind=%ld dateintmean = %lf dateintsum=%lf k2cpt=%lf k2=%lf\n",iind, dateintsum/k2cpt, dateintsum,k2cpt, k2); */
                                          }             }
                                          /*}*/             /*}*/
                                  } /* end m */           } /* end m */
                          } /* end bool */         } /* end bool */
                  } /* end iind = 1 to imx */       } /* end iind = 1 to imx */
        /* prop[s][age] is feeded for any initial and valid live state as well as         /* prop[s][age] is feeded for any initial and valid live state as well as
                                         freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */            freq[s1][s2][age] at single age of beginning the transition, for a combination j1 */
   
   
                  /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/       /*      fprintf(ficresp, "#Count between %.lf/%.lf/%.lf and %.lf/%.lf/%.lf\n",jprev1, mprev1,anprev1,jprev2, mprev2,anprev2);*/
                  pstamp(ficresp);       pstamp(ficresp);
                  if  (cptcovn>0) {       if  (cptcovn>0) {
                          fprintf(ficresp, "\n#********** Variable ");          fprintf(ficresp, "\n#********** Variable "); 
                          fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable ");          fprintf(ficresphtm, "\n<br/><br/><h3>********** Variable "); 
                          fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable ");          fprintf(ficresphtmfr, "\n<br/><br/><h3>********** Variable "); 
                          for (z1=1; z1<=cptcoveff; z1++){         for (z1=1; z1<=cptcoveff; z1++){
                                  fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);           fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
                                  fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);           fprintf(ficresphtm, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
                                  fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);           fprintf(ficresphtmfr, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
                          }         }
                          fprintf(ficresp, "**********\n#");         fprintf(ficresp, "**********\n#");
                          fprintf(ficresphtm, "**********</h3>\n");         fprintf(ficresphtm, "**********</h3>\n");
                          fprintf(ficresphtmfr, "**********</h3>\n");         fprintf(ficresphtmfr, "**********</h3>\n");
                          fprintf(ficlog, "\n#********** Variable ");          fprintf(ficlog, "\n#********** Variable "); 
                          for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficlog, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
                          fprintf(ficlog, "**********\n");         fprintf(ficlog, "**********\n");
                  }       }
                  fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");       fprintf(ficresphtm,"<table style=\"text-align:center; border: 1px solid\">");
                  for(i=1; i<=nlstate;i++) {       for(i=1; i<=nlstate;i++) {
                          fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);         fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
                          fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);         fprintf(ficresphtm, "<th>Age</th><th>Prev(%d)</th><th>N(%d)</th><th>N</th>",i,i);
                  }       }
                  fprintf(ficresp, "\n");       fprintf(ficresp, "\n");
                  fprintf(ficresphtm, "\n");       fprintf(ficresphtm, "\n");
               
                  /* Header of frequency table by age */       /* Header of frequency table by age */
                  fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");       fprintf(ficresphtmfr,"<table style=\"text-align:center; border: 1px solid\">");
                  fprintf(ficresphtmfr,"<th>Age</th> ");       fprintf(ficresphtmfr,"<th>Age</th> ");
                  for(jk=-1; jk <=nlstate+ndeath; jk++){       for(jk=-1; jk <=nlstate+ndeath; jk++){
                          for(m=-1; m <=nlstate+ndeath; m++){         for(m=-1; m <=nlstate+ndeath; m++){
                                  if(jk!=0 && m!=0)           if(jk!=0 && m!=0)
                                          fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);             fprintf(ficresphtmfr,"<th>%d%d</th> ",jk,m);
                          }         }
                  }       }
                  fprintf(ficresphtmfr, "\n");       fprintf(ficresphtmfr, "\n");
               
                  /* For each age */       /* For each age */
                  for(iage=iagemin; iage <= iagemax+3; iage++){       for(iage=iagemin; iage <= iagemax+3; iage++){
                          fprintf(ficresphtm,"<tr>");         fprintf(ficresphtm,"<tr>");
                          if(iage==iagemax+1){         if(iage==iagemax+1){
                                  fprintf(ficlog,"1");           fprintf(ficlog,"1");
                                  fprintf(ficresphtmfr,"<tr><th>0</th> ");           fprintf(ficresphtmfr,"<tr><th>0</th> ");
                          }else if(iage==iagemax+2){         }else if(iage==iagemax+2){
                                  fprintf(ficlog,"0");           fprintf(ficlog,"0");
                                  fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");           fprintf(ficresphtmfr,"<tr><th>Unknown</th> ");
                          }else if(iage==iagemax+3){         }else if(iage==iagemax+3){
                                  fprintf(ficlog,"Total");           fprintf(ficlog,"Total");
                                  fprintf(ficresphtmfr,"<tr><th>Total</th> ");           fprintf(ficresphtmfr,"<tr><th>Total</th> ");
                          }else{         }else{
                                  if(first==1){           if(first==1){
                                          first=0;             first=0;
                                          printf("See log file for details...\n");             printf("See log file for details...\n");
                                  }           }
                                  fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);           fprintf(ficresphtmfr,"<tr><th>%d</th> ",iage);
                                  fprintf(ficlog,"Age %d", iage);           fprintf(ficlog,"Age %d", iage);
                          }         }
                          for(jk=1; jk <=nlstate ; jk++){         for(jk=1; jk <=nlstate ; jk++){
                                  for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)           for(m=-1, pp[jk]=0; m <=nlstate+ndeath ; m++)
                                          pp[jk] += freq[jk][m][iage];              pp[jk] += freq[jk][m][iage]; 
                          }         }
                          for(jk=1; jk <=nlstate ; jk++){         for(jk=1; jk <=nlstate ; jk++){
                                  for(m=-1, pos=0; m <=0 ; m++)           for(m=-1, pos=0; m <=0 ; m++)
                                          pos += freq[jk][m][iage];             pos += freq[jk][m][iage];
                                  if(pp[jk]>=1.e-10){           if(pp[jk]>=1.e-10){
                                          if(first==1){             if(first==1){
                                                  printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);               printf(" %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
                                          }             }
                                          fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);             fprintf(ficlog," %d.=%.0f loss[%d]=%.1f%%",jk,pp[jk],jk,100*pos/pp[jk]);
                                  }else{           }else{
                                          if(first==1)             if(first==1)
                                                  printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);               printf(" %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
                                          fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);             fprintf(ficlog," %d.=%.0f loss[%d]=NaNQ%%",jk,pp[jk],jk);
                                  }           }
                          }         }
   
                          for(jk=1; jk <=nlstate ; jk++){          for(jk=1; jk <=nlstate ; jk++){ 
                                  /* posprop[jk]=0; */           /* posprop[jk]=0; */
                                  for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */           for(m=0, pp[jk]=0; m <=nlstate+ndeath; m++)/* Summing on all ages */
                                          pp[jk] += freq[jk][m][iage];             pp[jk] += freq[jk][m][iage];
                          }      /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */         }        /* pp[jk] is the total number of transitions starting from state jk and any ending status until this age */
   
                          for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){         for(jk=1,pos=0, pospropta=0.; jk <=nlstate ; jk++){
                                  pos += pp[jk]; /* pos is the total number of transitions until this age */           pos += pp[jk]; /* pos is the total number of transitions until this age */
                                  posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state           posprop[jk] += prop[jk][iage]; /* prop is the number of transitions from a live state
                                                                                                                                                                          from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */                                             from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
                                  pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state           pospropta += prop[jk][iage]; /* prop is the number of transitions from a live state
                                                                                                                                                                          from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */                                           from jk at age iage prop[s[m][iind]][(int)agev[m][iind]] += weight[iind] */
                          }         }
                          for(jk=1; jk <=nlstate ; jk++){         for(jk=1; jk <=nlstate ; jk++){
                                  if(pos>=1.e-5){           if(pos>=1.e-5){
                                          if(first==1)             if(first==1)
                                                  printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);               printf(" %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
                                          fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);             fprintf(ficlog," %d.=%.0f prev[%d]=%.1f%%",jk,pp[jk],jk,100*pp[jk]/pos);
                                  }else{           }else{
                                          if(first==1)             if(first==1)
                                                  printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);               printf(" %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
                                          fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);             fprintf(ficlog," %d.=%.0f prev[%d]=NaNQ%%",jk,pp[jk],jk);
                                  }           }
                                  if( iage <= iagemax){           if( iage <= iagemax){
                                          if(pos>=1.e-5){             if(pos>=1.e-5){
                                                  fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);               fprintf(ficresp," %d %.5f %.0f %.0f",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
                                                  fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);               fprintf(ficresphtm,"<th>%d</th><td>%.5f</td><td>%.0f</td><td>%.0f</td>",iage,prop[jk][iage]/pospropta, prop[jk][iage],pospropta);
                                                  /*probs[iage][jk][j1]= pp[jk]/pos;*/               /*probs[iage][jk][j1]= pp[jk]/pos;*/
                                                  /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/               /*printf("\niage=%d jk=%d j1=%d %.5f %.0f %.0f %f",iage,jk,j1,pp[jk]/pos, pp[jk],pos,probs[iage][jk][j1]);*/
                                          }             }
                                          else{             else{
                                                  fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);               fprintf(ficresp," %d NaNq %.0f %.0f",iage,prop[jk][iage],pospropta);
                                                  fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);               fprintf(ficresphtm,"<th>%d</th><td>NaNq</td><td>%.0f</td><td>%.0f</td>",iage, prop[jk][iage],pospropta);
                                          }             }
                                  }           }
                                  pospropt[jk] +=posprop[jk];           pospropt[jk] +=posprop[jk];
                          } /* end loop jk */         } /* end loop jk */
                          /* pospropt=0.; */         /* pospropt=0.; */
                          for(jk=-1; jk <=nlstate+ndeath; jk++){         for(jk=-1; jk <=nlstate+ndeath; jk++){
                                  for(m=-1; m <=nlstate+ndeath; m++){           for(m=-1; m <=nlstate+ndeath; m++){
                                          if(freq[jk][m][iage] !=0 ) { /* minimizing output */             if(freq[jk][m][iage] !=0 ) { /* minimizing output */
                                                  if(first==1){               if(first==1){
                                                          printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);                 printf(" %d%d=%.0f",jk,m,freq[jk][m][iage]);
                                                  }               }
                                                  fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);               fprintf(ficlog," %d%d=%.0f",jk,m,freq[jk][m][iage]);
                                          }             }
                                          if(jk!=0 && m!=0)             if(jk!=0 && m!=0)
                                                  fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);               fprintf(ficresphtmfr,"<td>%.0f</td> ",freq[jk][m][iage]);
                                  }           }
                          } /* end loop jk */         } /* end loop jk */
                          posproptt=0.;          posproptt=0.; 
                          for(jk=1; jk <=nlstate; jk++){         for(jk=1; jk <=nlstate; jk++){
                                  posproptt += pospropt[jk];           posproptt += pospropt[jk];
                          }         }
                          fprintf(ficresphtmfr,"</tr>\n ");         fprintf(ficresphtmfr,"</tr>\n ");
                          if(iage <= iagemax){         if(iage <= iagemax){
                                  fprintf(ficresp,"\n");           fprintf(ficresp,"\n");
                                  fprintf(ficresphtm,"</tr>\n");           fprintf(ficresphtm,"</tr>\n");
                          }         }
                          if(first==1)         if(first==1)
                                  printf("Others in log...\n");           printf("Others in log...\n");
                          fprintf(ficlog,"\n");         fprintf(ficlog,"\n");
                  } /* end loop age iage */       } /* end loop age iage */
                  fprintf(ficresphtm,"<tr><th>Tot</th>");       fprintf(ficresphtm,"<tr><th>Tot</th>");
                  for(jk=1; jk <=nlstate ; jk++){       for(jk=1; jk <=nlstate ; jk++){
                          if(posproptt < 1.e-5){         if(posproptt < 1.e-5){
                                  fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);             fprintf(ficresphtm,"<td>Nanq</td><td>%.0f</td><td>%.0f</td>",pospropt[jk],posproptt);  
                          }else{         }else{
                                  fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);              fprintf(ficresphtm,"<td>%.5f</td><td>%.0f</td><td>%.0f</td>",pospropt[jk]/posproptt,pospropt[jk],posproptt);   
                          }         }
                  }       }
                  fprintf(ficresphtm,"</tr>\n");       fprintf(ficresphtm,"</tr>\n");
                  fprintf(ficresphtm,"</table>\n");       fprintf(ficresphtm,"</table>\n");
                  fprintf(ficresphtmfr,"</table>\n");       fprintf(ficresphtmfr,"</table>\n");
                  if(posproptt < 1.e-5){       if(posproptt < 1.e-5){
                          fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);         fprintf(ficresphtm,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
                          fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);         fprintf(ficresphtmfr,"\n <p><b> This combination (%d) is not valid and no result will be produced</b></p>",j1);
                          fprintf(ficres,"\n  This combination (%d) is not valid and no result will be produced\n\n",j1);         fprintf(ficres,"\n  This combination (%d) is not valid and no result will be produced\n\n",j1);
                          invalidvarcomb[j1]=1;         invalidvarcomb[j1]=1;
                  }else{       }else{
                          fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);         fprintf(ficresphtm,"\n <p> This combination (%d) is valid and result will be produced.</p>",j1);
                          invalidvarcomb[j1]=0;         invalidvarcomb[j1]=0;
                  }       }
                  fprintf(ficresphtmfr,"</table>\n");       fprintf(ficresphtmfr,"</table>\n");
          } /* end selected combination of covariate j1 */     } /* end selected combination of covariate j1 */
          dateintmean=dateintsum/k2cpt;      dateintmean=dateintsum/k2cpt; 
                                     
          fclose(ficresp);     fclose(ficresp);
          fclose(ficresphtm);     fclose(ficresphtm);
          fclose(ficresphtmfr);     fclose(ficresphtmfr);
          free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE);     free_ma3x(freq,-5,nlstate+ndeath,-5,nlstate+ndeath, iagemin-AGEMARGE, iagemax+3+AGEMARGE);
          free_vector(pospropt,1,nlstate);     free_vector(pospropt,1,nlstate);
          free_vector(posprop,1,nlstate);     free_vector(posprop,1,nlstate);
          free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);     free_matrix(prop,1,nlstate,iagemin-AGEMARGE, iagemax+3+AGEMARGE);
          free_vector(pp,1,nlstate);     free_vector(pp,1,nlstate);
          /* End of Freq */     /* End of Freq */
  }   }
   
 /************ Prevalence ********************/  /************ Prevalence ********************/
Line 4104  void  concatwav(int wav[], int **dh, int Line 4191  void  concatwav(int wav[], int **dh, int
      and mw[mi+1][i]. dh depends on stepm.       and mw[mi+1][i]. dh depends on stepm.
      */       */
   
   int i, mi, m;    int i=0, mi=0, m=0, mli=0;
   /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;    /* int j, k=0,jk, ju, jl,jmin=1e+5, jmax=-1;
      double sum=0., jmean=0.;*/       double sum=0., jmean=0.;*/
   int first, firstwo, firsthree, firstfour;    int first=0, firstwo=0, firsthree=0, firstfour=0, firstfiv=0;
   int j, k=0,jk, ju, jl;    int j, k=0,jk, ju, jl;
   double sum=0.;    double sum=0.;
   first=0;  
   firstwo=0;  
   firsthree=0;  
   firstfour=0;  
   jmin=100000;    jmin=100000;
   jmax=-1;    jmax=-1;
   jmean=0.;    jmean=0.;
   
   /* Treating live states */
   for(i=1; i<=imx; i++){  /* For simple cases and if state is death */    for(i=1; i<=imx; i++){  /* For simple cases and if state is death */
     mi=0;      mi=0;  /* First valid wave */
                   mli=0; /* Last valid wave? */
     m=firstpass;      m=firstpass;
     while(s[m][i] <= nlstate){  /* a live state */      while(s[m][i] <= nlstate){  /* a live state */
       if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */                          if(m >firstpass && s[m][i]==s[m-1][i] && mint[m][i]==mint[m-1][i] && anint[m][i]==anint[m-1][i]){/* Two succesive identical information on wave m */
         mw[++mi][i]=m;                                  mli=m-1;/* mw[++mi][i]=m-1; */
       }                          }else if(s[m][i]>=1 || s[m][i]==-4 || s[m][i]==-5){ /* Since 0.98r4 if status=-2 vital status is really unknown, wave should be skipped */
       if(m >=lastpass){                                  mw[++mi][i]=m;
         if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){                                  mli=m;
           if(firsthree == 0){        } /* else might be a useless wave  -1 and mi is not incremented and mw[mi] not updated */
             printf("Information! Unknown health status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);        if(m < lastpass){ /* m < lastpass, standard case */
             firsthree=1;                                  m++; /* mi gives the "effective" current wave, m the current wave, go to next wave by incrementing m */
           }  
           fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);  
           mw[++mi][i]=m;  
         }  
         if(s[m][i]==-2){ /* Vital status is really unknown */  
           nbwarn++;  
           if((int)anint[m][i] == 9999){  /*  Has the vital status really been verified? */  
             printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);  
             fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);  
           }  
           break;  
         }  
         break;  
       }        }
       else                          else{ /* m >= lastpass, eventual special issue with warning */
         m++;  #ifdef UNKNOWNSTATUSNOTCONTRIBUTING
                                   break;
   #else
                                   if(s[m][i]==-1 && (int) andc[i] == 9999 && (int)anint[m][i] != 9999){
                                           if(firsthree == 0){
                                                   printf("Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
                                                   firsthree=1;
                                           }
                                           fprintf(ficlog,"Information! Unknown status for individual %ld line=%d occurred at last wave %d at known date %d/%d. Please, check if your unknown date of death %d/%d means a live state %d at wave %d. This case(%d)/wave(%d) contributes to the likelihood as pi. .\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], (int) moisdc[i], (int) andc[i], s[m][i], m, i, m);
                                           mw[++mi][i]=m;
                                           mli=m;
                                   }
                                   if(s[m][i]==-2){ /* Vital status is really unknown */
                                           nbwarn++;
                                           if((int)anint[m][i] == 9999){  /*  Has the vital status really been verified? */
                                                   printf("Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
                                                   fprintf(ficlog,"Warning! Vital status for individual %ld (line=%d) at last wave %d interviewed at date %d/%d is unknown %d. Please, check if the vital status and the date of death %d/%d are really unknown. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], (int) moisdc[i], (int) andc[i], i, m);
                                           }
                                           break;
                                   }
                                   break;
   #endif
                           }/* End m >= lastpass */
     }/* end while */      }/* end while */
       
           /* mi is the last effective wave, m is lastpass, mw[j][i] gives the # of j-th effective wave for individual i */
     /* After last pass */      /* After last pass */
   /* Treating death states */
     if (s[m][i] > nlstate){  /* In a death state */      if (s[m][i] > nlstate){  /* In a death state */
                           /* if( mint[m][i]==mdc[m][i] && anint[m][i]==andc[m][i]){ /\* same date of death and date of interview *\/ */
                           /* } */
       mi++;     /* Death is another wave */        mi++;     /* Death is another wave */
       /* if(mi==0)  never been interviewed correctly before death */        /* if(mi==0)  never been interviewed correctly before death */
          /* Only death is a correct wave */                          /* Only death is a correct wave */
       mw[mi][i]=m;        mw[mi][i]=m;
     }else if ((int) andc[i] != 9999) { /* Status is either death or negative. A death occured after lastpass, we can't take it into account because of potential bias */      }
   #ifndef DISPATCHINGKNOWNDEATHAFTERLASTWAVE
                   else if ((int) andc[i] != 9999) { /* Status is negative. A death occured after lastpass, we can't take it into account because of potential bias */
       /* m++; */        /* m++; */
       /* mi++; */        /* mi++; */
       /* s[m][i]=nlstate+1;  /\* We are setting the status to the last of non live state *\/ */        /* s[m][i]=nlstate+1;  /\* We are setting the status to the last of non live state *\/ */
       /* mw[mi][i]=m; */        /* mw[mi][i]=m; */
       nberr++;  
       if ((int)anint[m][i]!= 9999) { /* date of last interview is known */        if ((int)anint[m][i]!= 9999) { /* date of last interview is known */
         if(firstwo==0){                                  if((andc[i]+moisdc[i]/12.) <=(anint[m][i]+mint[m][i]/12.)){ /* death occured before last wave and status should have been death instead of -1 */
           printf("Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );                                          nbwarn++;
           firstwo=1;                                          if(firstfiv==0){
         }                                                  printf("Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m );
         fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );                                                  firstfiv=1;
                                           }else{
                                                   fprintf(ficlog,"Warning! Death for individual %ld line=%d occurred at %d/%d before last wave %d interviewed at %d/%d and should have been coded as death instead of '%d'. This case (%d)/wave (%d) is contributing to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], s[m][i], i,m );
                                           }
                                   }else{ /* Death occured afer last wave potential bias */
                                           nberr++;
                                           if(firstwo==0){
                                                   printf("Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
                                                   firstwo=1;
                                           }
                                           fprintf(ficlog,"Error! Death for individual %ld line=%d occurred at %d/%d after last wave %d interviewed at %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
                                   }
       }else{ /* end date of interview is known */        }else{ /* end date of interview is known */
         /* death is known but not confirmed by death status at any wave */                                  /* death is known but not confirmed by death status at any wave */
         if(firstfour==0){                                  if(firstfour==0){
           printf("Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );                                          printf("Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\nOthers in log file only\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
           firstfour=1;                                          firstfour=1;
         }                                  }
         fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );                                  fprintf(ficlog,"Error! Death for individual %ld line=%d  occurred %d/%d but not confirmed by any death status for any wave, including last wave %d at unknown date %d/%d. Potential bias if other individuals are still alive at this date but ignored. This case (%d)/wave (%d) is skipped, no contribution to likelihood.\n",num[i],i,(int) moisdc[i], (int) andc[i], lastpass,(int)mint[m][i],(int)anint[m][i], i,m );
       }        }
     }      } /* end if date of death is known */
     wav[i]=mi;  #endif
       wav[i]=mi; /* mi should be the last effective wave (or mli) */
       /* wav[i]=mw[mi][i]; */
     if(mi==0){      if(mi==0){
       nbwarn++;        nbwarn++;
       if(first==0){        if(first==0){
         printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);                                  printf("Warning! No valid information for individual %ld line=%d (skipped) and may be others, see log file\n",num[i],i);
         first=1;                                  first=1;
       }        }
       if(first==1){        if(first==1){
         fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);                                  fprintf(ficlog,"Warning! No valid information for individual %ld line=%d (skipped)\n",num[i],i);
       }        }
     } /* end mi==0 */      } /* end mi==0 */
   } /* End individuals */    } /* End individuals */
   /* wav and mw are no more changed */    /* wav and mw are no more changed */
           
       
   for(i=1; i<=imx; i++){    for(i=1; i<=imx; i++){
     for(mi=1; mi<wav[i];mi++){      for(mi=1; mi<wav[i];mi++){
       if (stepm <=0)        if (stepm <=0)
         dh[mi][i]=1;                                  dh[mi][i]=1;
       else{        else{
         if (s[mw[mi+1][i]][i] > nlstate) { /* A death */                                  if (s[mw[mi+1][i]][i] > nlstate) { /* A death */
           if (agedc[i] < 2*AGESUP) {                                          if (agedc[i] < 2*AGESUP) {
             j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12);                                                   j= rint(agedc[i]*12-agev[mw[mi][i]][i]*12); 
             if(j==0) j=1;  /* Survives at least one month after exam */                                                  if(j==0) j=1;  /* Survives at least one month after exam */
             else if(j<0){                                                  else if(j<0){
               nberr++;                                                          nberr++;
               printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                                                          printf("Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               j=1; /* Temporary Dangerous patch */                                                          j=1; /* Temporary Dangerous patch */
               printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);                                                          printf("   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
               fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                                                          fprintf(ficlog,"Error! Negative delay (%d to death) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
               fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);                                                          fprintf(ficlog,"   We assumed that the date of interview was correct (and not the date of death) and postponed the death %d month(s) (one stepm) after the interview. You MUST fix the contradiction between dates.\n",stepm);
             }                                                  }
             k=k+1;                                                  k=k+1;
             if (j >= jmax){                                                  if (j >= jmax){
               jmax=j;                                                          jmax=j;
               ijmax=i;                                                          ijmax=i;
             }                                                  }
             if (j <= jmin){                                                  if (j <= jmin){
               jmin=j;                                                          jmin=j;
               ijmin=i;                                                          ijmin=i;
             }                                                  }
             sum=sum+j;                                                  sum=sum+j;
             /*if (j<0) printf("j=%d num=%d \n",j,i);*/                                                  /*if (j<0) printf("j=%d num=%d \n",j,i);*/
             /*    printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/                                                  /*        printf("%d %d %d %d\n", s[mw[mi][i]][i] ,s[mw[mi+1][i]][i],j,i);*/
           }                                          }
         }                                  }
         else{                                  else{
           j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));                                          j= rint( (agev[mw[mi+1][i]][i]*12 - agev[mw[mi][i]][i]*12));
 /*        if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */  /*        if (j<0) printf("%d %lf %lf %d %d %d\n", i,agev[mw[mi+1][i]][i], agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]); */
                                           
           k=k+1;                                          k=k+1;
           if (j >= jmax) {                                          if (j >= jmax) {
             jmax=j;                                                  jmax=j;
             ijmax=i;                                                  ijmax=i;
           }                                          }
           else if (j <= jmin){                                          else if (j <= jmin){
             jmin=j;                                                  jmin=j;
             ijmin=i;                                                  ijmin=i;
           }                                          }
           /*        if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */                                          /*          if (j<10) printf("j=%d jmin=%d num=%d ",j,jmin,i); */
           /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/                                          /*printf("%d %lf %d %d %d\n", i,agev[mw[mi][i]][i],j,s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);*/
           if(j<0){                                          if(j<0){
             nberr++;                                                  nberr++;
             printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                                                  printf("Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
             fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);                                                  fprintf(ficlog,"Error! Negative delay (%d) between waves %d and %d of individual %ld at line %d who is aged %.1f with statuses from %d to %d\n ",j,mw[mi][i],mw[mi+1][i],num[i], i,agev[mw[mi][i]][i],s[mw[mi][i]][i] ,s[mw[mi+1][i]][i]);
           }                                          }
           sum=sum+j;                                          sum=sum+j;
         }                                  }
         jk= j/stepm;                                  jk= j/stepm;
         jl= j -jk*stepm;                                  jl= j -jk*stepm;
         ju= j -(jk+1)*stepm;                                  ju= j -(jk+1)*stepm;
         if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */                                  if(mle <=1){ /* only if we use a the linear-interpoloation pseudo-likelihood */
           if(jl==0){                                          if(jl==0){
             dh[mi][i]=jk;                                                  dh[mi][i]=jk;
             bh[mi][i]=0;                                                  bh[mi][i]=0;
           }else{ /* We want a negative bias in order to only have interpolation ie                                          }else{ /* We want a negative bias in order to only have interpolation ie
                   * to avoid the price of an extra matrix product in likelihood */                                                                          * to avoid the price of an extra matrix product in likelihood */
             dh[mi][i]=jk+1;                                                  dh[mi][i]=jk+1;
             bh[mi][i]=ju;                                                  bh[mi][i]=ju;
           }                                          }
         }else{                                  }else{
           if(jl <= -ju){                                          if(jl <= -ju){
             dh[mi][i]=jk;                                                  dh[mi][i]=jk;
             bh[mi][i]=jl;       /* bias is positive if real duration                                                  bh[mi][i]=jl;   /* bias is positive if real duration
                                  * is higher than the multiple of stepm and negative otherwise.                                                                                                           * is higher than the multiple of stepm and negative otherwise.
                                  */                                                                                                           */
           }                                          }
           else{                                          else{
             dh[mi][i]=jk+1;                                                  dh[mi][i]=jk+1;
             bh[mi][i]=ju;                                                  bh[mi][i]=ju;
           }                                          }
           if(dh[mi][i]==0){                                          if(dh[mi][i]==0){
             dh[mi][i]=1; /* At least one step */                                                  dh[mi][i]=1; /* At least one step */
             bh[mi][i]=ju; /* At least one step */                                                  bh[mi][i]=ju; /* At least one step */
             /*  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/                                                  /*  printf(" bh=%d ju=%d jl=%d dh=%d jk=%d stepm=%d %d\n",bh[mi][i],ju,jl,dh[mi][i],jk,stepm,i);*/
           }                                          }
         } /* end if mle */                                  } /* end if mle */
       }        }
     } /* end wave */      } /* end wave */
   }    }
Line 4345  void  concatwav(int wav[], int **dh, int Line 4458  void  concatwav(int wav[], int **dh, int
                                                                                                                                 undefined. Usually 3: -1, 0 and 1. */                                                                                                                                  undefined. Usually 3: -1, 0 and 1. */
       }        }
       /* In fact  ncodemax[j]=2 (dichotom. variables only) but it could be more for        /* In fact  ncodemax[j]=2 (dichotom. variables only) but it could be more for
                                  historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */                           * historical reasons: 3 if coded 1, 2, 3 and 4 and Ndum[2]=0 */
     } /* Ndum[-1] number of undefined modalities */      } /* Ndum[-1] number of undefined modalities */
                                   
     /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */      /* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
Line 4366  void  concatwav(int wav[], int **dh, int Line 4479  void  concatwav(int wav[], int **dh, int
         if (Ndum[i] == 0) { /* If nobody responded to this modality k */          if (Ndum[i] == 0) { /* If nobody responded to this modality k */
                                 break;                                  break;
                         }                          }
         ij++;                          ij++;
         nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/                          nbcode[Tvar[j]][ij]=i;  /* stores the original value of modality i in an array nbcode, ij modality from 1 to last non-nul modality.*/
         cptcode = ij; /* New max modality for covar j */                          cptcode = ij; /* New max modality for covar j */
     } /* end of loop on modality i=-1 to 1 or more */      } /* end of loop on modality i=-1 to 1 or more */
                         
     /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */      /*   for (k=0; k<= cptcode; k++) { /\* k=-1 ? k=0 to 1 *\//\* Could be 1 to 4 *\//\* cptcode=modmaxcovj *\/ */
     /*  /\*recode from 0 *\/ */      /*  /\*recode from 0 *\/ */
     /*                               k is a modality. If we have model=V1+V1*sex  */      /*                               k is a modality. If we have model=V1+V1*sex  */
Line 4392  void  concatwav(int wav[], int **dh, int Line 4505  void  concatwav(int wav[], int **dh, int
                 /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/                   /* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/ 
                 ij=Tvar[i]; /* Tvar might be -1 if status was unknown */                   ij=Tvar[i]; /* Tvar might be -1 if status was unknown */ 
                 Ndum[ij]++; /* Might be supersed V1 + V1*age */                  Ndum[ij]++; /* Might be supersed V1 + V1*age */
         }           } /* V4+V3+V5, Ndum[1]@5={0, 0, 1, 1, 1} */
                   
         ij=0;          ij=0;
         for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */          for (i=0; i<=  maxncov-1; i++) { /* modmaxcovj is unknown here. Only Ndum[2(V2),3(age*V3), 5(V3*V2) 6(V1*V4) */
Line 4407  void  concatwav(int wav[], int **dh, int Line 4520  void  concatwav(int wav[], int **dh, int
         }          }
         /* ij--; */          /* ij--; */
         /* cptcoveff=ij; /\*Number of total covariates*\/ */          /* cptcoveff=ij; /\*Number of total covariates*\/ */
         *cptcov=ij; /*Number of total covariates*/          *cptcov=ij; /*Number of total real effective covariates: effective
                                                                    * because they can be excluded from the model and real
                                                            * if in the model but excluded because missing values*/
 }  }
   
   
Line 5259  To be simple, these graphs help to under Line 5373  To be simple, these graphs help to under
    tj = (int) pow(2,cptcoveff);     tj = (int) pow(2,cptcoveff);
    if (cptcovn<1) {tj=1;ncodemax[1]=1;}     if (cptcovn<1) {tj=1;ncodemax[1]=1;}
    j1=0;     j1=0;
    for(j1=1; j1<=tj;j1++){  /* For each valid combination of covariates */     for(j1=1; j1<=tj;j1++){  /* For each valid combination of covariates or only once*/
      if  (cptcovn>0) {       if  (cptcovn>0) {
        fprintf(ficresprob, "\n#********** Variable ");          fprintf(ficresprob, "\n#********** Variable "); 
        for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);         for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresprob, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtabm(j1,z1)]);
Line 5700  true period expectancies (those weighted Line 5814  true period expectancies (those weighted
 }  }
   
 /******************* Gnuplot file **************/  /******************* Gnuplot file **************/
  void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){  void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double fage , int prevfcast, int backcast, char pathc[], double p[]){
   
   char dirfileres[132],optfileres[132];    char dirfileres[132],optfileres[132];
   char gplotcondition[132];    char gplotcondition[132];
Line 5963  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6077  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* Survival functions (period) from state i in state j by final state j */    /* Survival functions (period) from state i in state j by final state j */
   for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */    for (k1=1; k1<= m ; k1 ++) { /* For each covariate if any */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each inital state  */
                           
       fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);        fprintf(ficgp,"\n#\n#\n# Survival functions in state j and all livestates from state i by final state j: 'lij' files, cov=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */        for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */          lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
Line 5975  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6089  plot [%.f:%.f]  ", ageminpar, agemaxpar)
       }        }
       fprintf(ficgp,"\n#\n");        fprintf(ficgp,"\n#\n");
       if(invalidvarcomb[k1]){        if(invalidvarcomb[k1]){
         fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);                                   fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
         continue;                                  continue;
       }        }
                                                   
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"LIJT_"),cpt,k1);
Line 5986  unset log y\n Line 6100  unset log y\n
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3;        k=3;
       for (j=1; j<= nlstate ; j ++){ /* Lived in state j */        for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
         if(j==1)                                  if(j==1)
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));                                          fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
         else                                  else
           fprintf(ficgp,", '' ");                                          fprintf(ficgp,", '' ");
         l=(nlstate+ndeath)*(cpt-1) +j;                                  l=(nlstate+ndeath)*(cpt-1) +j;
         fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);                                  fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):($%d",k1,k+l);
         /* for (i=2; i<= nlstate+ndeath ; i ++) */                                  /* for (i=2; i<= nlstate+ndeath ; i ++) */
         /*   fprintf(ficgp,"+$%d",k+l+i-1); */                                  /*   fprintf(ficgp,"+$%d",k+l+i-1); */
         fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);                                  fprintf(ficgp,") t \"l(%d,%d)\" w l",cpt,j);
       } /* nlstate */        } /* nlstate */
       fprintf(ficgp,", '' ");        fprintf(ficgp,", '' ");
       fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);        fprintf(ficgp," u (($1==%d && (floor($2)%%5 == 0)) ? ($3):1/0):(",k1);
       for (j=1; j<= nlstate ; j ++){ /* Lived in state j */        for (j=1; j<= nlstate ; j ++){ /* Lived in state j */
         l=(nlstate+ndeath)*(cpt-1) +j;                                  l=(nlstate+ndeath)*(cpt-1) +j;
         if(j < nlstate)                                  if(j < nlstate)
           fprintf(ficgp,"$%d +",k+l);                                          fprintf(ficgp,"$%d +",k+l);
         else                                  else
           fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);                                          fprintf(ficgp,"$%d) t\"l(%d,.)\" w l",k+l,cpt);
       }        }
       fprintf(ficgp,"\nset out\n");        fprintf(ficgp,"\nset out\n");
     } /* end cpt state*/       } /* end cpt state*/ 
Line 6013  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6127  plot [%.f:%.f]  ", ageminpar, agemaxpar)
   /* CV preval stable (period) for each covariate */    /* CV preval stable (period) for each covariate */
   for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */    for (k1=1; k1<= m ; k1 ++) { /* For each covariate combination (1 to m=2**k), if any covariate is present */
     for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */      for (cpt=1; cpt<=nlstate ; cpt ++) { /* For each life state */
                           
       fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);        fprintf(ficgp,"\n#\n#\n#CV preval stable (period): 'pij' files, covariatecombination#=%d state=%d",k1, cpt);
       for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */        for (k=1; k<=cptcoveff; k++){    /* For each covariate and each value */
         lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */          lv= decodtabm(k1,k,cptcoveff); /* Should be the covariate number corresponding to k1 combination */
Line 6025  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6139  plot [%.f:%.f]  ", ageminpar, agemaxpar)
       }        }
       fprintf(ficgp,"\n#\n");        fprintf(ficgp,"\n#\n");
       if(invalidvarcomb[k1]){        if(invalidvarcomb[k1]){
         fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);                                   fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
         continue;                                  continue;
       }        }
                           
       fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);        fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"P_"),cpt,k1);
       fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\        fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
 set ter svg size 640, 480\n\  set ter svg size 640, 480\n                                                                                                                                                                              \
 unset log y\n\  unset log y\n                                                                                                                                                                                                                                    \
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
       k=3; /* Offset */        k=3; /* Offset */
       for (i=1; i<= nlstate ; i ++){        for (i=1; i<= nlstate ; i ++){
         if(i==1)                                  if(i==1)
           fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));                                          fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJ_"));
         else                                  else
           fprintf(ficgp,", '' ");                                          fprintf(ficgp,", '' ");
         l=(nlstate+ndeath)*(i-1)+1;                                  l=(nlstate+ndeath)*(i-1)+1;
         fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);                                  fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l);
         for (j=2; j<= nlstate ; j ++)                                  for (j=2; j<= nlstate ; j ++)
           fprintf(ficgp,"+$%d",k+l+j-1);                                          fprintf(ficgp,"+$%d",k+l+j-1);
         fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);                                  fprintf(ficgp,")) t \"prev(%d,%d)\" w l",i,cpt);
       } /* nlstate */        } /* nlstate */
       fprintf(ficgp,"\nset out\n");        fprintf(ficgp,"\nset out\n");
     } /* end cpt state*/       } /* end cpt state*/ 
   } /* end covariate */      } /* end covariate */  
           
           
 /* 7eme */  /* 7eme */
   if(backcast == 1){    if(backcast == 1){
     /* CV back preval stable (period) for each covariate */      /* CV back preval stable (period) for each covariate */
Line 6062  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6176  plot [%.f:%.f]  ", ageminpar, agemaxpar)
           /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */            /* decodtabm(1,1,4) = 1 because h=1  k= (1) 1  1  1 */
           /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */            /* decodtabm(1,2,4) = 1 because h=1  k=  1 (1) 1  1 */
           /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */            /* decodtabm(13,3,4)= 2 because h=13 k=  1  1 (2) 2 */
           vlv= nbcode[Tvaraff[k]][lv];                                          vlv= nbcode[Tvaraff[k]][lv];
           fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);                                          fprintf(ficgp," V%d=%d ",Tvaraff[k],vlv);
         }                                  }
         fprintf(ficgp,"\n#\n");                                  fprintf(ficgp,"\n#\n");
         if(invalidvarcomb[k1]){                                  if(invalidvarcomb[k1]){
           fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1);                                           fprintf(ficgp,"#Combination (%d) ignored because no cases \n",k1); 
           continue;                                          continue;
         }                                  }
                                                                   
         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);                                  fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PB_"),cpt,k1);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\                                  fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Probability\" \n\
 set ter svg size 640, 480\n                                                                                                                                                                                     \  set ter svg size 640, 480\n                                                                                                                                                                                     \
 unset log y\n                                                                                                                                                                                                                                           \  unset log y\n                                                                                                                                                                                                                                           \
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
         k=3; /* Offset */                                  k=3; /* Offset */
         for (i=1; i<= nlstate ; i ++){                                  for (i=1; i<= nlstate ; i ++){
           if(i==1)                                          if(i==1)
             fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));                                                  fprintf(ficgp,"\"%s\"",subdirf2(fileresu,"PIJB_"));
           else                                          else
             fprintf(ficgp,", '' ");                                                  fprintf(ficgp,", '' ");
           /* l=(nlstate+ndeath)*(i-1)+1; */                                          /* l=(nlstate+ndeath)*(i-1)+1; */
           l=(nlstate+ndeath)*(cpt-1)+1;                                          l=(nlstate+ndeath)*(cpt-1)+1;
           /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */                                          /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l); /\* a vérifier *\/ */
           /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */                                          /* fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d/($%d",k1,k+l+(cpt-1),k+l+(cpt-1)+i-1); /\* a vérifier *\/ */
           fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */                                          fprintf(ficgp," u ($1==%d ? ($3):1/0):($%d",k1,k+l+(cpt-1)+i-1); /* a vérifier */
           /* for (j=2; j<= nlstate ; j ++) */                                          /* for (j=2; j<= nlstate ; j ++) */
           /*    fprintf(ficgp,"+$%d",k+l+j-1); */                                          /*      fprintf(ficgp,"+$%d",k+l+j-1); */
           /*    /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */                                          /*      /\* fprintf(ficgp,"+$%d",k+l+j-1); *\/ */
           fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);                                          fprintf(ficgp,") t \"bprev(%d,%d)\" w l",i,cpt);
         } /* nlstate */                                  } /* nlstate */
         fprintf(ficgp,"\nset out\n");                                  fprintf(ficgp,"\nset out\n");
       } /* end cpt state*/         } /* end cpt state*/ 
     } /* end covariate */        } /* end covariate */  
   } /* End if backcast */    } /* End if backcast */
Line 6118  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6232  plot [%.f:%.f]  ", ageminpar, agemaxpar)
           continue;            continue;
         }          }
                                                                   
         fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");                                  fprintf(ficgp,"# hpijx=probability over h years, hp.jx is weighted by observed prev\n ");
         fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);                                  fprintf(ficgp,"\nset out \"%s_%d-%d.svg\" \n",subdirf2(optionfilefiname,"PROJ_"),cpt,k1);
         fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\                                  fprintf(ficgp,"set xlabel \"Age\" \nset ylabel \"Prevalence\" \n\
 set ter svg size 640, 480\n     \  set ter svg size 640, 480\n                                                                                                                                                                                     \
 unset log y\n   \  unset log y\n                                                                                                                                                                                                                                           \
 plot [%.f:%.f]  ", ageminpar, agemaxpar);  plot [%.f:%.f]  ", ageminpar, agemaxpar);
         for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */          for (i=1; i<= nlstate+1 ; i ++){  /* nlstate +1 p11 p21 p.1 */
           /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/            /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
Line 6196  plot [%.f:%.f]  ", ageminpar, agemaxpar) Line 6310  plot [%.f:%.f]  ", ageminpar, agemaxpar)
     fprintf(ficgp,"# initial state %d\n",i);      fprintf(ficgp,"# initial state %d\n",i);
     for(k=1; k <=(nlstate+ndeath); k++){      for(k=1; k <=(nlstate+ndeath); k++){
       if (k != i) {        if (k != i) {
         fprintf(ficgp,"#   current state %d\n",k);                                  fprintf(ficgp,"#   current state %d\n",k);
         for(j=1; j <=ncovmodel; j++){                                  for(j=1; j <=ncovmodel; j++){
           fprintf(ficgp,"p%d=%f; ",jk,p[jk]);                                          fprintf(ficgp,"p%d=%f; ",jk,p[jk]);
           jk++;                                           jk++; 
         }                                  }
         fprintf(ficgp,"\n");                                  fprintf(ficgp,"\n");
       }        }
     }      }
   }    }
   fprintf(ficgp,"##############\n#\n");    fprintf(ficgp,"##############\n#\n");
           
   /*goto avoid;*/    /*goto avoid;*/
   fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");    fprintf(ficgp,"\n##############\n#Graphics of probabilities or incidences\n#############\n");
   fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");    fprintf(ficgp,"# logi(p12/p11)=a12+b12*age+c12age*age+d12*V1+e12*V1*age\n");
Line 7196  int readdata(char datafile[], int firsto Line 7310  int readdata(char datafile[], int firsto
     }      }
     trimbb(linetmp,line); /* Trims multiple blanks in line */      trimbb(linetmp,line); /* Trims multiple blanks in line */
     strcpy(line, linetmp);      strcpy(line, linetmp);
         
 /* Loops on waves */      /* Loops on waves */
     for (j=maxwav;j>=1;j--){      for (j=maxwav;j>=1;j--){
       cutv(stra, strb, line, ' ');         cutv(stra, strb, line, ' '); 
       if(strb[0]=='.') { /* Missing status */        if(strb[0]=='.') { /* Missing status */
Line 7212  int readdata(char datafile[], int firsto Line 7326  int readdata(char datafile[], int firsto
           return 1;            return 1;
         }          }
       }        }
        
       s[j][i]=lval;        s[j][i]=lval;
         
         /* Date of Interview */
       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{
         printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);                                  printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);
         fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);                                  fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d.  Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
         return 1;                                  return 1;
       }        }
       anint[j][i]= (double) year;         anint[j][i]= (double) year; 
       mint[j][i]= (double)month;         mint[j][i]= (double)month; 
Line 7240  int readdata(char datafile[], int firsto Line 7356  int readdata(char datafile[], int firsto
       year=9999;        year=9999;
     }else{      }else{
       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);        printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
         fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);                          fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of death (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
         return 1;                          return 1;
     }      }
     andc[i]=(double) year;       andc[i]=(double) year; 
     moisdc[i]=(double) month;       moisdc[i]=(double) month; 
Line 7257  int readdata(char datafile[], int firsto Line 7373  int readdata(char datafile[], int firsto
     }else{      }else{
       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);        printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);
       fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);        fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .).  Exiting.\n",strb, linei,i,line);fflush(ficlog);
         return 1;                          return 1;
     }      }
     if (year==9999) {      if (year==9999) {
       printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given.  Exiting.\n",strb, linei,i,line);        printf("Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given.  Exiting.\n",strb, linei,i,line);
       fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);        fprintf(ficlog,"Error reading data around '%s' at line number %d for individual %d, '%s'\nShould be a date of birth (mm/yyyy) but at least the year of birth should be given. Exiting.\n",strb, linei,i,line);fflush(ficlog);
         return 1;                          return 1;
   
     }      }
     annais[i]=(double)(year);      annais[i]=(double)(year);
Line 7286  int readdata(char datafile[], int firsto Line 7402  int readdata(char datafile[], int firsto
     for (j=ncovcol;j>=1;j--){      for (j=ncovcol;j>=1;j--){
       cutv(stra, strb,line,' ');         cutv(stra, strb,line,' '); 
       if(strb[0]=='.') { /* Missing covariate value */        if(strb[0]=='.') { /* Missing covariate value */
         lval=-1;                                  lval=-1;
       }else{        }else{
         errno=0;                                  errno=0;
         lval=strtol(strb,&endptr,10);                                   lval=strtol(strb,&endptr,10); 
         if( strb[0]=='\0' || (*endptr != '\0')){                                  if( strb[0]=='\0' || (*endptr != '\0')){
           printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);                                          printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);
           fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);                                          fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative).  Exiting.\n",lval, linei,i, line);fflush(ficlog);
           return 1;                                          return 1;
         }                                  }
       }        }
       if(lval <-1 || lval >1){        if(lval <-1 || lval >1){
         printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \                                  printf("Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \   Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \   for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
  For example, for multinomial values like 1, 2 and 3,\n \   For example, for multinomial values like 1, 2 and 3,\n \
Line 7306  int readdata(char datafile[], int firsto Line 7422  int readdata(char datafile[], int firsto
  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \   and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
  output of IMaCh is often meaningless.\n \   output of IMaCh is often meaningless.\n \
  Exiting.\n",lval,linei, i,line,j);   Exiting.\n",lval,linei, i,line,j);
         fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \                                  fprintf(ficlog,"Error reading data around '%ld' at line number %d for individual %d, '%s'\n \
  Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \   Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
  for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \   for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
  For example, for multinomial values like 1, 2 and 3,\n \   For example, for multinomial values like 1, 2 and 3,\n \
Line 7315  int readdata(char datafile[], int firsto Line 7431  int readdata(char datafile[], int firsto
  and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \   and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
  output of IMaCh is often meaningless.\n \   output of IMaCh is often meaningless.\n \
  Exiting.\n",lval,linei, i,line,j);fflush(ficlog);   Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
         return 1;                                  return 1;
       }        }
       covar[j][i]=(double)(lval);        covar[j][i]=(double)(lval);
       strcpy(line,stra);        strcpy(line,stra);
Line 7339  int readdata(char datafile[], int firsto Line 7455  int readdata(char datafile[], int firsto
     
   return (0);    return (0);
   /* endread: */    /* endread: */
     printf("Exiting readdata: ");          printf("Exiting readdata: ");
     fclose(fic);          fclose(fic);
     return (1);          return (1);
   
   
   
 }  }
   
 void removespace(char *str) {  void removespace(char *str) {
   char *p1 = str, *p2 = str;    char *p1 = str, *p2 = str;
   do    do
Line 7354  void removespace(char *str) { Line 7468  void removespace(char *str) {
   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)
    * Model  V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age   /**< This routine decode the model and returns:
    * - nagesqr = 1 if age*age in the model, otherwise 0.          * Model  V1+V2+V3+V8+V7*V8+V5*V6+V8*age+V3*age+age*age
    * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age          * - nagesqr = 1 if age*age in the model, otherwise 0.
    * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age          * - cptcovt total number of covariates of the model nbocc(+)+1 = 8 excepting constant and age and age*age
    * - cptcovage number of covariates with age*products =2          * - cptcovn or number of covariates k of the models excluding age*products =6 and age*age
    * - cptcovs number of simple covariates          * - cptcovage number of covariates with age*products =2
    * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10          * - cptcovs number of simple covariates
    *     which is a new column after the 9 (ncovcol) variables.           * - Tvar[k] is the id of the kth covariate Tvar[1]@12 {1, 2, 3, 8, 10, 11, 8, 3, 7, 8, 5, 6}, thus Tvar[5=V7*V8]=10
    * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual          *     which is a new column after the 9 (ncovcol) variables. 
    * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage          * - if k is a product Vn*Vm covar[k][i] is filled with correct values for each individual
    *    Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.          * - Tprod[l] gives the kth covariates of the product Vn*Vm l=1 to cptcovprod-cptcovage
    * - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .          *    Tprod[1]@2 {5, 6}: position of first product V7*V8 is 5, and second V5*V6 is 6.
  */          * - Tvard[k]  p Tvard[1][1]@4 {7, 8, 5, 6} for V7*V8 and V5*V6 .
           */
 {  {
   int i, j, k, ks;    int i, j, k, ks;
   int  j1, k1, k2;    int  j1, k1, k2;
Line 7392  int decodemodel ( char model[], int last Line 7507  int decodemodel ( char model[], int last
     if ((strpt=strstr(model,"age*age")) !=0){      if ((strpt=strstr(model,"age*age")) !=0){
       printf(" strpt=%s, model=%s\n",strpt, model);        printf(" strpt=%s, model=%s\n",strpt, model);
       if(strpt != model){        if(strpt != model){
       printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \                                  printf("Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \   'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model);   corresponding column of parameters.\n",model);
       fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \                                  fprintf(ficlog,"Error in model: 'model=%s'; 'age*age' should in first place before other covariates\n \
  'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \   'model=1+age+age*age+V1.' or 'model=1+age+age*age+V1+V1*age.', please swap as well as \n \
  corresponding column of parameters.\n",model); fflush(ficlog);   corresponding column of parameters.\n",model); fflush(ficlog);
       return 1;                                  return 1;
     }                          }
   
       nagesqr=1;        nagesqr=1;
       if (strstr(model,"+age*age") !=0)        if (strstr(model,"+age*age") !=0)
         substrchaine(modelsav, model, "+age*age");                                  substrchaine(modelsav, model, "+age*age");
       else if (strstr(model,"age*age+") !=0)        else if (strstr(model,"age*age+") !=0)
         substrchaine(modelsav, model, "age*age+");                                  substrchaine(modelsav, model, "age*age+");
       else         else 
         substrchaine(modelsav, model, "age*age");                                  substrchaine(modelsav, model, "age*age");
     }else      }else
       nagesqr=0;        nagesqr=0;
     if (strlen(modelsav) >1){      if (strlen(modelsav) >1){
       j=nbocc(modelsav,'+'); /**< j=Number of '+' */        j=nbocc(modelsav,'+'); /**< j=Number of '+' */
       j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */        j1=nbocc(modelsav,'*'); /**< j1=Number of '*' */
       cptcovs=j+1-j1; /**<  Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =2  */        cptcovs=j+1-j1; /**<  Number of simple covariates V1+V1*age+V3 +V3*V4+age*age=> V1 + V3 =5-3=2  */
       cptcovt= j+1; /* Number of total covariates in the model, not including        cptcovt= j+1; /* Number of total covariates in the model, not including
                    * cst, age and age*age                                                                                    * cst, age and age*age 
                    * V1+V1*age+ V3 + V3*V4+age*age=> 4*/                                                                                   * V1+V1*age+ V3 + V3*V4+age*age=> 3+1=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                            * but the covariates which are products must be treated 
                   * separately: ncovn=4- 2=2 (V1+V3). */                           * 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  */
   
Line 7431  int decodemodel ( char model[], int last Line 7546  int decodemodel ( char model[], int last
        *   k=  1    2      3       4     5       6      7        8         *   k=  1    2      3       4     5       6      7        8
        *  cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8         *  cptcovn number of covariates (not including constant and age ) = # of + plus 1 = 7+1=8
        *  covar[k,i], value of kth covariate if not including age for individual i:         *  covar[k,i], value of kth covariate if not including age for individual i:
        *       covar[1][i]= (V2), covar[4][i]=(V3), covar[8][i]=(V8)         *       covar[1][i]= (V1), covar[4][i]=(V4), covar[8][i]=(V8)
        *  Tvar[k] # of the kth covariate:  Tvar[1]=2  Tvar[4]=3 Tvar[8]=8         *  Tvar[k] # of the kth covariate:  Tvar[1]=2  Tvar[2]=1 Tvar[4]=3 Tvar[8]=8
        *       if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and          *       if multiplied by age: V3*age Tvar[3=V3*age]=3 (V3) Tvar[7]=8 and 
        *  Tage[++cptcovage]=k         *  Tage[++cptcovage]=k
        *       if products, new covar are created after ncovcol with k1         *       if products, new covar are created after ncovcol with k1
Line 7476  int decodemodel ( char model[], int last Line 7591  int decodemodel ( char model[], int last
         Tvar[k]=0;          Tvar[k]=0;
       cptcovage=0;        cptcovage=0;
       for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */        for(k=1; k<=cptcovt;k++){ /* Loop on total covariates of the model */
         cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+'                                   cutl(stra,strb,modelsav,'+'); /* keeps in strb after the first '+' 
                                          modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */                                                                                                                                                                    modelsav==V2+V1+V4+V3*age strb=V3*age stra=V2+V1+V4 */ 
         if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */                                  if (nbocc(modelsav,'+')==0) strcpy(strb,modelsav); /* and analyzes it */
         /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/                                  /*      printf("i=%d a=%s b=%s sav=%s\n",i, stra,strb,modelsav);*/
         /*scanf("%d",i);*/                                  /*scanf("%d",i);*/
         if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */                                  if (strchr(strb,'*')) {  /**< Model includes a product V2+V1+V4+V3*age strb=V3*age */
           cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */                                          cutl(strc,strd,strb,'*'); /**< strd*strc  Vm*Vn: strb=V3*age(input) strc=age strd=V3 ; V3*V2 strc=V2, strd=V3 */
           if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */                                          if (strcmp(strc,"age")==0) { /**< Model includes age: Vn*age */
             /* covar is not filled and then is empty */                                                  /* covar is not filled and then is empty */
             cptcovprod--;                                                  cptcovprod--;
             cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */                                                  cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
             Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */                                                  Tvar[k]=atoi(stre);  /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2; V1+V1*age Tvar[2]=1 */
             cptcovage++; /* Sums the number of covariates which include age as a product */                                                  cptcovage++; /* Sums the number of covariates which include age as a product */
             Tage[cptcovage]=k;  /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */                                                  Tage[cptcovage]=k;  /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
             /*printf("stre=%s ", stre);*/                                                  /*printf("stre=%s ", stre);*/
           } else if (strcmp(strd,"age")==0) { /* or age*Vn */                                          } else if (strcmp(strd,"age")==0) { /* or age*Vn */
             cptcovprod--;                                                  cptcovprod--;
             cutl(stre,strb,strc,'V');                                                  cutl(stre,strb,strc,'V');
             Tvar[k]=atoi(stre);                                                  Tvar[k]=atoi(stre);
             cptcovage++;                                                  cptcovage++;
             Tage[cptcovage]=k;                                                  Tage[cptcovage]=k;
           } else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/                                          } else {  /* Age is not in the model product V2+V1+V1*V4+V3*age+V3*V2  strb=V3*V2*/
             /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */                                                  /* loops on k1=1 (V3*V2) and k1=2 V4*V3 */
             cptcovn++;                                                  cptcovn++;
             cptcovprodnoage++;k1++;                                                  cptcovprodnoage++;k1++;
             cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/                                                  cutl(stre,strb,strc,'V'); /* strc= Vn, stre is n; strb=V3*V2 stre=3 strc=*/
             Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but                                                  Tvar[k]=ncovcol+k1; /* For model-covariate k tells which data-covariate to use but
                                    because this model-covariate is a construction we invent a new column                                                                                                                                           because this model-covariate is a construction we invent a new column
                                    ncovcol + k1                                                                                                                                           ncovcol + k1
                                    If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2                                                                                                                                           If already ncovcol=4 and model=V2+V1+V1*V4+age*V3+V3*V2
                                    Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */                                                                                                                                           Tvar[3=V1*V4]=4+1 Tvar[5=V3*V2]=4 + 2= 6, etc */
             cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */                                                  cutl(strc,strb,strd,'V'); /* strd was Vm, strc is m */
             Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */                                                  Tprod[k1]=k;  /* Tprod[1]=3(=V1*V4) for V2+V1+V1*V4+age*V3+V3*V2  */
             Tvard[k1][1] =atoi(strc); /* m 1 for V1*/                                                  Tvard[k1][1] =atoi(strc); /* m 1 for V1*/
             Tvard[k1][2] =atoi(stre); /* n 4 for V4*/                                                  Tvard[k1][2] =atoi(stre); /* n 4 for V4*/
             k2=k2+2;                                                  k2=k2+2;
             Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */                                                  Tvar[cptcovt+k2]=Tvard[k1][1]; /* Tvar[(cptcovt=4+k2=1)=5]= 1 (V1) */
             Tvar[cptcovt+k2+1]=Tvard[k1][2];  /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */                                                  Tvar[cptcovt+k2+1]=Tvard[k1][2];  /* Tvar[(cptcovt=4+(k2=1)+1)=6]= 4 (V4) */
             for (i=1; i<=lastobs;i++){                                                  for (i=1; i<=lastobs;i++){
               /* Computes the new covariate which is a product of                                                          /* Computes the new covariate which is a product of
                  covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */                                                                   covar[n][i]* covar[m][i] and stores it at ncovol+k1 May not be defined */
               covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];                                                          covar[ncovcol+k1][i]=covar[atoi(stre)][i]*covar[atoi(strc)][i];
             }                                                  }
           } /* End age is not in the model */                                          } /* End age is not in the model */
         } /* End if model includes a product */                                  } /* End if model includes a product */
         else { /* no more sum */                                  else { /* no more sum */
           /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/                                          /*printf("d=%s c=%s b=%s\n", strd,strc,strb);*/
           /*  scanf("%d",i);*/                                          /*  scanf("%d",i);*/
           cutl(strd,strc,strb,'V');                                          cutl(strd,strc,strb,'V');
           ks++; /**< Number of simple covariates */                                          ks++; /**< Number of simple covariates */
           cptcovn++;                                          cptcovn++;
           Tvar[k]=atoi(strd);                                          Tvar[k]=atoi(strd);
         }                                  }
         strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */                                   strcpy(modelsav,stra);  /* modelsav=V2+V1+V4 stra=V2+V1+V4 */ 
         /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);                                  /*printf("a=%s b=%s sav=%s\n", stra,strb,modelsav);
           scanf("%d",i);*/                                          scanf("%d",i);*/
       } /* end of loop + on total covariates */        } /* end of loop + on total covariates */
     } /* end if strlen(modelsave == 0) age*age might exist */      } /* end if strlen(modelsave == 0) age*age might exist */
   } /* end if strlen(model == 0) */    } /* end if strlen(model == 0) */
Line 7540  int decodemodel ( char model[], int last Line 7655  int decodemodel ( char model[], int last
     If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/      If model=V1+V1*age then Tvar[1]=1 Tvar[2]=1 cptcovage=1 Tage[1]=2 cptcovprod=0*/
   
   /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);    /* printf("tvar1=%d tvar2=%d tvar3=%d cptcovage=%d Tage=%d",Tvar[1],Tvar[2],Tvar[3],cptcovage,Tage[1]);
   printf("cptcovprod=%d ", cptcovprod);                   printf("cptcovprod=%d ", cptcovprod);
   fprintf(ficlog,"cptcovprod=%d ", cptcovprod);                   fprintf(ficlog,"cptcovprod=%d ", cptcovprod);
   
   scanf("%d ",i);*/                   scanf("%d ",i);*/
   /* Dispatching in quantitative and time varying covariates */
   
   
   return (0); /* with covar[new additional covariate if product] and Tage if age */     return (0); /* with covar[new additional covariate if product] and Tage if age */ 
   /*endread:*/    /*endread:*/
     printf("Exiting decodemodel: ");          printf("Exiting decodemodel: ");
     return (1);          return (1);
 }  }
   
 int 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 )
Line 8742  Please run with mle=-1 to get a correct Line 8858  Please run with mle=-1 to get a correct
       printf("Problem writing new parameter file: %s\n", rfileres);goto end;        printf("Problem writing new parameter file: %s\n", rfileres);goto end;
       fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;        fprintf(ficlog,"Problem writing new parameter file: %s\n", rfileres);goto end;
     }      }
     fprintf(ficres,"#%s\n",version);      fprintf(ficres,"#IMaCh version %s\n",version);
   }    /* End of mle != -3 */    }    /* End of mle != -3 */
       
   /*  Main data    /*  Main data
Line 8810  Please run with mle=-1 to get a correct Line 8926  Please run with mle=-1 to get a correct
 /* Main decodemodel */  /* Main decodemodel */
   
   
   if(decodemodel(model, lastobs) == 1)    if(decodemodel(model, lastobs) == 1) /* In order to get Tvar[k] V4+V3+V5 p Tvar[1]@3  = {4, 3, 5}*/
     goto end;      goto end;
   
   if((double)(lastobs-imx)/(double)imx > 1.10){    if((double)(lastobs-imx)/(double)imx > 1.10){
Line 9067  Interval (in months) between two waves: Line 9183  Interval (in months) between two waves:
     ageexmed=vector(1,n);      ageexmed=vector(1,n);
     agecens=vector(1,n);      agecens=vector(1,n);
     dcwave=ivector(1,n);      dcwave=ivector(1,n);
                    
     for (i=1; i<=imx; i++){      for (i=1; i<=imx; i++){
       dcwave[i]=-1;        dcwave[i]=-1;
       for (m=firstpass; m<=lastpass; m++)        for (m=firstpass; m<=lastpass; m++)
Line 9567  Please run with mle=-1 to get a correct Line 9683  Please run with mle=-1 to get a correct
     ungetc(c,ficpar);      ungetc(c,ficpar);
           
     fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);      fscanf(ficpar,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",&backcast,&jback1,&mback1,&anback1,&jback2,&mback2,&anback2,&mobilavproj);
     fprintf(ficparo,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);      fprintf(ficparo,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
     fprintf(ficlog,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);      fprintf(ficlog,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
     fprintf(ficres,"backcast=%d starting-back-date=%lf/%lf/%lf final-back-date=%lf/%lf/%lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);      fprintf(ficres,"backcast=%d starting-back-date=%.lf/%.lf/%.lf final-back-date=%.lf/%.lf/%.lf mobil_average=%d\n",backcast,jback1,mback1,anback1,jback2,mback2,anback2,mobilavproj);
     /* day and month of proj2 are not used but only year anproj2.*/      /* day and month of proj2 are not used but only year anproj2.*/
           
           
Removed from v.1.1  
changed lines
  Added in v.1.2

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