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

Diff for /imach/src/imach.c between versions 1.167 and 1.181

version 1.167, 2014/12/22 13:50:56 version 1.181, 2015/02/11 23:22:24
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.181  2015/02/11 23:22:24  brouard
     Summary: Comments on Powell added
   
     Author:
   
     Revision 1.180  2015/02/11 17:33:45  brouard
     Summary: Finishing move from main to function (hpijx and prevalence_limit)
   
     Revision 1.179  2015/01/04 09:57:06  brouard
     Summary: back to OS/X
   
     Revision 1.178  2015/01/04 09:35:48  brouard
     *** empty log message ***
   
     Revision 1.177  2015/01/03 18:40:56  brouard
     Summary: Still testing ilc32 on OSX
   
     Revision 1.176  2015/01/03 16:45:04  brouard
     *** empty log message ***
   
     Revision 1.175  2015/01/03 16:33:42  brouard
     *** empty log message ***
   
     Revision 1.174  2015/01/03 16:15:49  brouard
     Summary: Still in cross-compilation
   
     Revision 1.173  2015/01/03 12:06:26  brouard
     Summary: trying to detect cross-compilation
   
     Revision 1.172  2014/12/27 12:07:47  brouard
     Summary: Back from Visual Studio and Intel, options for compiling for Windows XP
   
     Revision 1.171  2014/12/23 13:26:59  brouard
     Summary: Back from Visual C
   
     Still problem with utsname.h on Windows
   
     Revision 1.170  2014/12/23 11:17:12  brouard
     Summary: Cleaning some \%% back to %%
   
     The escape was mandatory for a specific compiler (which one?), but too many warnings.
   
     Revision 1.169  2014/12/22 23:08:31  brouard
     Summary: 0.98p
   
     Outputs some informations on compiler used, OS etc. Testing on different platforms.
   
     Revision 1.168  2014/12/22 15:17:42  brouard
     Summary: update
   
   Revision 1.167  2014/12/22 13:50:56  brouard    Revision 1.167  2014/12/22 13:50:56  brouard
   Summary: Testing uname and compiler version and if compiled 32 or 64    Summary: Testing uname and compiler version and if compiled 32 or 64
   
Line 511 Line 561
 */  */
   
 #define POWELL /* Instead of NLOPT */  #define POWELL /* Instead of NLOPT */
   #define POWELLDIRECT /* Directest to decide new direction instead of Powell test */
   
 #include <math.h>  #include <math.h>
 #include <stdio.h>  #include <stdio.h>
Line 519 Line 570
   
 #ifdef _WIN32  #ifdef _WIN32
 #include <io.h>  #include <io.h>
   #include <windows.h>
   #include <tchar.h>
 #else  #else
 #include <unistd.h>  #include <unistd.h>
 #endif  #endif
   
 #include <limits.h>  #include <limits.h>
 #include <sys/types.h>  #include <sys/types.h>
 #include <sys/utsname.h>  
   #if defined(__GNUC__)
   #include <sys/utsname.h> /* Doesn't work on Windows */
   #endif
   
 #include <sys/stat.h>  #include <sys/stat.h>
 #include <errno.h>  #include <errno.h>
 /* extern int errno; */  /* extern int errno; */
Line 590  typedef struct { Line 647  typedef struct {
 /* $Id$ */  /* $Id$ */
 /* $State$ */  /* $State$ */
   
 char version[]="Imach version 0.99, September 2014,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121)";  char version[]="Imach version 0.98p, Février 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
 char fullversion[]="$Revision$ $Date$";   char fullversion[]="$Revision$ $Date$"; 
 char strstart[80];  char strstart[80];
 char optionfilext[10], optionfilefiname[FILENAMELENGTH];  char optionfilext[10], optionfilefiname[FILENAMELENGTH];
Line 1353  void powell(double p[], double **xi, int Line 1410  void powell(double p[], double **xi, int
               double (*func)(double []));                 double (*func)(double [])); 
   int i,ibig,j;     int i,ibig,j; 
   double del,t,*pt,*ptt,*xit;    double del,t,*pt,*ptt,*xit;
     double directest;
   double fp,fptt;    double fp,fptt;
   double *xits;    double *xits;
   int niterf, itmp;    int niterf, itmp;
Line 1413  void powell(double p[], double **xi, int Line 1471  void powell(double p[], double **xi, int
       printf("%d",i);fflush(stdout);        printf("%d",i);fflush(stdout);
       fprintf(ficlog,"%d",i);fflush(ficlog);        fprintf(ficlog,"%d",i);fflush(ficlog);
       linmin(p,xit,n,fret,func);         linmin(p,xit,n,fret,func); 
       if (fabs(fptt-(*fret)) > del) {         if (fabs(fptt-(*fret)) > del) { /* We are keeping the max gain on each of the n directions 
                                          because that direction will be replaced unless the gain del is small
                                         in comparison with the 'probable' gain, mu^2, with the last average direction.
                                         Unless the n directions are conjugate some gain in the determinant may be obtained
                                         with the new direction.
                                         */
         del=fabs(fptt-(*fret));           del=fabs(fptt-(*fret)); 
         ibig=i;           ibig=i; 
       }         } 
Line 1465  void powell(double p[], double **xi, int Line 1528  void powell(double p[], double **xi, int
       return;         return; 
     }       } 
     if (*iter == ITMAX) nrerror("powell exceeding maximum iterations.");       if (*iter == ITMAX) nrerror("powell exceeding maximum iterations."); 
     for (j=1;j<=n;j++) { /* Computes an extrapolated point */      for (j=1;j<=n;j++) { /* Computes the extrapolated point P_0 + 2 (P_n-P_0) */
       ptt[j]=2.0*p[j]-pt[j];         ptt[j]=2.0*p[j]-pt[j]; 
       xit[j]=p[j]-pt[j];         xit[j]=p[j]-pt[j]; 
       pt[j]=p[j];         pt[j]=p[j]; 
     }       } 
     fptt=(*func)(ptt);       fptt=(*func)(ptt); /* f_3 */
     if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */      if (fptt < fp) { /* If extrapolated point is better, decide if we keep that new direction or not */
       /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */        /* (x1 f1=fp), (x2 f2=*fret), (x3 f3=fptt), (xm fm) */
       /* From x1 (P0) distance of x2 is at h and x3 is 2h */        /* From x1 (P0) distance of x2 is at h and x3 is 2h */
       /* Let f"(x2) be the 2nd derivative equal everywhere.  */        /* Let f"(x2) be the 2nd derivative equal everywhere.  */
       /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */        /* Then the parabolic through (x1,f1), (x2,f2) and (x3,f3) */
       /* will reach at f3 = fm + h^2/2 f"m  ; f" = (f1 -2f2 +f3 ) / h**2 */        /* will reach at f3 = fm + h^2/2 f"m  ; f" = (f1 -2f2 +f3 ) / h**2 */
       /* f1-f3 = delta(2h) = 2 h**2 f'' = 2(f1- 2f2 +f3) */        /* Conditional for using this new direction is that mu^2 = (f1-2f2+f3)^2 /2 < del */
       /* Thus we compare delta(2h) with observed f1-f3 */  
       /* or best gain on one ancient line 'del' with total  */  
       /* gain f1-f2 = f1 - f2 - 'del' with del  */  
       /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */        /* t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt); */
   
       t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del);        t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del);
       t= t- del*SQR(fp-fptt);        t= t- del*SQR(fp-fptt);
       printf("t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);        directest = SQR(fp-2.0*(*fret)+fptt) - 2.0 * del; /* If del was big enough we change it for a new direction */
       fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t);  
 #ifdef DEBUG  #ifdef DEBUG
         printf("t1= %.12lf, t2= %.12lf, t=%.12lf  directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
         fprintf(ficlog,"t1= %.12lf, t2= %.12lf, t=%.12lf directest=%.12lf\n", 2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del),del*SQR(fp-fptt),t,directest);
       printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),        printf("t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
              (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));               (fp-(*fret)-del)*(fp-(*fret)-del),(fp-fptt)*(fp-fptt));
       fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),        fprintf(ficlog,"t3= %.12lf, t4= %.12lf, t3*= %.12lf, t4*= %.12lf\n",SQR(fp-(*fret)-del),SQR(fp-fptt),
Line 1495  void powell(double p[], double **xi, int Line 1556  void powell(double p[], double **xi, int
       printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);        printf("tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
       fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);        fprintf(ficlog, "tt= %.12lf, t=%.12lf\n",2.0*(fp-2.0*(*fret)+fptt)*(fp-(*fret)-del)*(fp-(*fret)-del)-del*(fp-fptt)*(fp-fptt),t);
 #endif  #endif
       if (t < 0.0) { /* Then we use it for last direction */  #ifdef POWELLDIRECT
         linmin(p,xit,n,fret,func); /* computes mean on the extrapolated direction.*/        if (directest < 0.0) { /* Then we use it for new direction */
   #else
         if (t < 0.0) { /* Then we use it for new direction */
   #endif
           linmin(p,xit,n,fret,func); /* computes minimum on the extrapolated direction.*/
         for (j=1;j<=n;j++) {           for (j=1;j<=n;j++) { 
           xi[j][ibig]=xi[j][n]; /* Replace the direction with biggest decrease by n */            xi[j][ibig]=xi[j][n]; /* Replace direction with biggest decrease by last direction n */
           xi[j][n]=xit[j];      /* and nth direction by the extrapolated */            xi[j][n]=xit[j];      /* and this nth direction by the by the average p_0 p_n */
         }          }
         printf("Gaining to use average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);          printf("Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
         fprintf(ficlog,"Gaining to use average direction of P0 P%d instead of biggest increase direction :\n",n,ibig);          fprintf(ficlog,"Gaining to use new average direction of P0 P%d instead of biggest increase direction %d :\n",n,ibig);
   
 #ifdef DEBUG  #ifdef DEBUG
         printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);          printf("Direction changed  last moved %d in place of ibig=%d, new last is the average:\n",n,ibig);
Line 1525  double **prevalim(double **prlim, int nl Line 1590  double **prevalim(double **prlim, int nl
 {  {
   /* Computes the prevalence limit in each live state at age x by left multiplying the unit    /* Computes the prevalence limit in each live state at age x by left multiplying the unit
      matrix by transitions matrix until convergence is reached */       matrix by transitions matrix until convergence is reached */
     
   int i, ii,j,k;    int i, ii,j,k;
   double min, max, maxmin, maxmax,sumnew=0.;    double min, max, maxmin, maxmax,sumnew=0.;
   /* double **matprod2(); */ /* test */    /* double **matprod2(); */ /* test */
   double **out, cov[NCOVMAX+1], **pmij();    double **out, cov[NCOVMAX+1], **pmij();
   double **newm;    double **newm;
   double agefin, delaymax=50 ; /* Max number of years to converge */    double agefin, delaymax=50 ; /* Max number of years to converge */
     
   for (ii=1;ii<=nlstate+ndeath;ii++)    for (ii=1;ii<=nlstate+ndeath;ii++)
     for (j=1;j<=nlstate+ndeath;j++){      for (j=1;j<=nlstate+ndeath;j++){
       oldm[ii][j]=(ii==j ? 1.0 : 0.0);        oldm[ii][j]=(ii==j ? 1.0 : 0.0);
     }      }
     
    cov[1]=1.;    cov[1]=1.;
      
  /* Even if hstepm = 1, at least one multiplication by the unit matrix */    /* Even if hstepm = 1, at least one multiplication by the unit matrix */
   for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){    for(agefin=age-stepm/YEARM; agefin>=age-delaymax; agefin=agefin-stepm/YEARM){
     newm=savm;      newm=savm;
     /* Covariates have to be included here again */      /* Covariates have to be included here again */
Line 1577  double **prevalim(double **prlim, int nl Line 1642  double **prevalim(double **prlim, int nl
       }        }
       maxmin=max-min;        maxmin=max-min;
       maxmax=FMAX(maxmax,maxmin);        maxmax=FMAX(maxmax,maxmin);
     }      } /* j loop */
     if(maxmax < ftolpl){      if(maxmax < ftolpl){
       return prlim;        return prlim;
     }      }
   }    } /* age loop */
     return prlim; /* should not reach here */
 }  }
   
 /*************** transition probabilities ***************/   /*************** transition probabilities ***************/ 
Line 2231  void mlikeli(FILE *ficres,double p[], in Line 2297  void mlikeli(FILE *ficres,double p[], in
 #endif  #endif
   free_matrix(xi,1,npar,1,npar);    free_matrix(xi,1,npar,1,npar);
   fclose(ficrespow);    fclose(ficrespow);
   printf("\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));    printf("#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
   fprintf(ficlog,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));    fprintf(ficlog,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
   fprintf(ficres,"\n#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));    fprintf(ficres,"#Number of iterations & function calls = %d & %d, -2 Log likelihood = %.12f\n",iter, countcallfunc,func(p));
   
 }  }
   
Line 2530  void  freqsummary(char fileres[], int ia Line 2596  void  freqsummary(char fileres[], int ia
   
   first=1;    first=1;
   
   /* for(k1=1; k1<=j ; k1++){   /* Loop on covariates */    /* for(k1=1; k1<=j ; k1++){ */  /* Loop on covariates */
   /*  for(i1=1; i1<=ncodemax[k1];i1++){ /* Now it is 2 */    /*  for(i1=1; i1<=ncodemax[k1];i1++){ */ /* Now it is 2 */
   /*    j1++;    /*    j1++; */
 */  
   for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){    for (j1 = 1; j1 <= (int) pow(2,cptcoveff); j1++){
       /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);        /*printf("cptcoveff=%d Tvaraff=%d", cptcoveff,Tvaraff[1]);
         scanf("%d", i);*/          scanf("%d", i);*/
Line 2910  void tricode(int *Tvar, int **nbcode, in Line 2975  void tricode(int *Tvar, int **nbcode, in
 {  {
   /**< Uses cptcovn+2*cptcovprod as the number of covariates */    /**< Uses cptcovn+2*cptcovprod as the number of covariates */
   /*      Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1     /*      Tvar[i]=atoi(stre);  find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 
   /* Boring subroutine which should only output nbcode[Tvar[j]][k]     * Boring subroutine which should only output nbcode[Tvar[j]][k]
    * Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2)     * Tvar[5] in V2+V1+V3*age+V2*V4 is 2 (V2)
   /* nbcode[Tvar[j]][1]=      * nbcode[Tvar[j]][1]= 
   */    */
   
   int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;    int ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
Line 3340  void varevsij(char optionfilefiname[], d Line 3405  void varevsij(char optionfilefiname[], d
   /* Variance of health expectancies */    /* Variance of health expectancies */
   /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/    /*  double **prevalim(double **prlim, int nlstate, double *xp, double age, double **oldm, double ** savm,double ftolpl);*/
   /* double **newm;*/    /* double **newm;*/
     /* int movingaverage(double ***probs, double bage,double fage, double ***mobaverage, int mobilav)*/
     
     int movingaverage();
   double **dnewm,**doldm;    double **dnewm,**doldm;
   double **dnewmp,**doldmp;    double **dnewmp,**doldmp;
   int i, j, nhstepm, hstepm, h, nstepm ;    int i, j, nhstepm, hstepm, h, nstepm ;
Line 3617  void varevsij(char optionfilefiname[], d Line 3685  void varevsij(char optionfilefiname[], d
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */  /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
 /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */  /*   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
   fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));    fprintf(ficgp,"\n plot \"%s\"  u 1:($3) not w l lt 1 ",subdirf(fileresprobmorprev));
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95\%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3+1.96*$4)) t \"95%% interval\" w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));    fprintf(ficgp,"\n replot \"%s\"  u 1:(($3-1.96*$4)) not w l lt 2 ",subdirf(fileresprobmorprev));
   fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));    fprintf(fichtm,"\n<br> File (multiple files are possible if covariates are present): <A href=\"%s\">%s</a>\n",subdirf(fileresprobmorprev),subdirf(fileresprobmorprev));
   fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);    fprintf(fichtm,"\n<br> Probability is computed over estepm=%d months. <br> <img src=\"%s%s.png\"> <br>\n", estepm,subdirf3(optionfilefiname,"varmuptjgr",digitp),digit);
Line 4038  To be simple, these graphs help to under Line 4106  To be simple, these graphs help to under
           } /* k12 */            } /* k12 */
         } /*l1 */          } /*l1 */
       }/* k1 */        }/* k1 */
       /* } /* loop covariates */        /* } */ /* loop covariates */
   }    }
   free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);    free_ma3x(varpij,1,nlstate,1,nlstate+ndeath,(int) bage, (int)fage);
   free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);    free_matrix(mu,1,(nlstate+ndeath)*(nlstate+ndeath),(int) bage, (int)fage);
Line 4215  void printinggnuplot(char fileres[], cha Line 4283  void printinggnuplot(char fileres[], cha
      fprintf(ficgp,"set xlabel \"Age\" \n\       fprintf(ficgp,"set xlabel \"Age\" \n\
 set ylabel \"Probability\" \n\  set ylabel \"Probability\" \n\
 set ter png small size 320, 240\n\  set ter png small size 320, 240\n\
 plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"\%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);  plot [%.f:%.f] \"%s\" every :::%d::%d u 1:2 \"%%lf",ageminpar,fage,subdirf2(fileres,"vpl"),k1-1,k1-1);
   
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else        fprintf(ficgp," \%%*lf (\%%*lf)");         else        fprintf(ficgp," %%*lf (%%*lf)");
      }       }
      fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);       fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," %%*lf (%%*lf)");
      }        } 
      fprintf(ficgp,"\" t\"95\%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);        fprintf(ficgp,"\" t\"95%% CI\" w l lt 1,\"%s\" every :::%d::%d u 1:($2-1.96*$3) \"%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1); 
      for (i=1; i<= nlstate ; i ++) {       for (i=1; i<= nlstate ; i ++) {
        if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");         if (i==cpt) fprintf(ficgp," %%lf (%%lf)");
        else fprintf(ficgp," \%%*lf (\%%*lf)");         else fprintf(ficgp," %%*lf (%%*lf)");
      }         }  
      fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));       fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" every :::%d::%d u 1:($%d) t\"Observed prevalence \" w l lt 2",subdirf2(fileres,"p"),k1-1,k1-1,2+4*(cpt-1));
    }     }
Line 4242  plot [%.f:%.f] \"%s\" every :::%d::%d u Line 4310  plot [%.f:%.f] \"%s\" every :::%d::%d u
           
     for (i=1; i<= nlstate+1 ; i ++) {      for (i=1; i<= nlstate+1 ; i ++) {
       k=2*i;        k=2*i;
       fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);        fprintf(ficgp,"\"%s\" every :::%d::%d u 1:2 \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
       for (j=1; j<= nlstate+1 ; j ++) {        for (j=1; j<= nlstate+1 ; j ++) {
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");          if (j==i) fprintf(ficgp," %%lf (%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");          else fprintf(ficgp," %%*lf (%%*lf)");
       }           }   
       if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");        if (i== 1) fprintf(ficgp,"\" t\"TLE\" w l ,");
       else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);        else fprintf(ficgp,"\" t\"LE in state (%d)\" w l ,",i-1);
       fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);        fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2-$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
       for (j=1; j<= nlstate+1 ; j ++) {        for (j=1; j<= nlstate+1 ; j ++) {
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");          if (j==i) fprintf(ficgp," %%lf (%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");          else fprintf(ficgp," %%*lf (%%*lf)");
       }           }   
       fprintf(ficgp,"\" t\"\" w l lt 0,");        fprintf(ficgp,"\" t\"\" w l lt 0,");
       fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"\%%lf",subdirf2(fileres,"t"),k1-1,k1-1);        fprintf(ficgp,"\"%s\" every :::%d::%d u 1:($2+$3*2) \"%%lf",subdirf2(fileres,"t"),k1-1,k1-1);
       for (j=1; j<= nlstate+1 ; j ++) {        for (j=1; j<= nlstate+1 ; j ++) {
         if (j==i) fprintf(ficgp," \%%lf (\%%lf)");          if (j==i) fprintf(ficgp," %%lf (%%lf)");
         else fprintf(ficgp," \%%*lf (\%%*lf)");          else fprintf(ficgp," %%*lf (%%*lf)");
       }           }   
       if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");        if (i== (nlstate+1)) fprintf(ficgp,"\" t\"\" w l lt 0");
       else fprintf(ficgp,"\" t\"\" w l lt 0,");        else fprintf(ficgp,"\" t\"\" w l lt 0,");
Line 4426  int movingaverage(double ***probs, doubl Line 4494  int movingaverage(double ***probs, doubl
   
   
 /************** Forecasting ******************/  /************** Forecasting ******************/
 prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){  void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anproj2, double p[], int cptcoveff){
   /* proj1, year, month, day of starting projection     /* proj1, year, month, day of starting projection 
      agemin, agemax range of age       agemin, agemax range of age
      dateprev1 dateprev2 range of dates during which prevalence is computed       dateprev1 dateprev2 range of dates during which prevalence is computed
Line 4549  prevforecast(char fileres[], double anpr Line 4617  prevforecast(char fileres[], double anpr
 }  }
   
 /************** Forecasting *****not tested NB*************/  /************** Forecasting *****not tested NB*************/
 populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){  void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){
       
   int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;    int cpt, stepsize, hstepm, nhstepm, j,k,c, cptcod, i,h;
   int *popage;    int *popage;
Line 5025  int readdata(char datafile[], int firsto Line 5093  int readdata(char datafile[], int firsto
               
       strcpy(line,stra);        strcpy(line,stra);
       cutv(stra, strb,line,' ');        cutv(stra, strb,line,' ');
       if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){        if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
       }        }
       else  if(iout=sscanf(strb,"%s.",dummy) != 0){        else  if( (iout=sscanf(strb,"%s.",dummy)) != 0){
         month=99;          month=99;
         year=9999;          year=9999;
       }else{        }else{
Line 5041  int readdata(char datafile[], int firsto Line 5109  int readdata(char datafile[], int firsto
     } /* ENd Waves */      } /* ENd Waves */
           
     cutv(stra, strb,line,' ');       cutv(stra, strb,line,' '); 
     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){      if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
     }      }
     else  if(iout=sscanf(strb,"%s.",dummy) != 0){      else  if( (iout=sscanf(strb,"%s.",dummy)) != 0){
       month=99;        month=99;
       year=9999;        year=9999;
     }else{      }else{
Line 5056  int readdata(char datafile[], int firsto Line 5124  int readdata(char datafile[], int firsto
     strcpy(line,stra);      strcpy(line,stra);
           
     cutv(stra, strb,line,' ');       cutv(stra, strb,line,' '); 
     if(iout=sscanf(strb,"%d/%d",&month, &year) != 0){      if( (iout=sscanf(strb,"%d/%d",&month, &year)) != 0){
     }      }
     else  if(iout=sscanf(strb,"%s.", dummy) != 0){      else  if( (iout=sscanf(strb,"%s.", dummy)) != 0){
       month=99;        month=99;
       year=9999;        year=9999;
     }else{      }else{
Line 5156  void removespace(char *str) { Line 5224  void removespace(char *str) {
   do    do
     while (*p2 == ' ')      while (*p2 == ' ')
       p2++;        p2++;
   while (*p1++ = *p2++);    while (*p1++ == *p2++);
 }  }
   
 int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:  int decodemodel ( char model[], int lastobs) /**< This routine decode the model and returns:
Line 5186  int decodemodel ( char model[], int last Line 5254  int decodemodel ( char model[], int last
     cptcovs=j+1-j1; /**<  Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2  */      cptcovs=j+1-j1; /**<  Number of simple covariates V1+V2*age+V3 +V3*V4=> V1 + V3 =2  */
     cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/      cptcovt= j+1; /* Number of total covariates in the model V1 + V2*age+ V3 + V3*V4=> 4*/
                   /* including age products which are counted in cptcovage.                    /* including age products which are counted in cptcovage.
                   /* but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */                    * but the covariates which are products must be treated separately: ncovn=4- 2=2 (V1+V3). */
     cptcovprod=j1; /**< Number of products  V1*V2 +v3*age = 2 */      cptcovprod=j1; /**< Number of products  V1*V2 +v3*age = 2 */
     cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1  */      cptcovprodnoage=0; /**< Number of covariate products without age: V3*V4 =1  */
     strcpy(modelsav,model);       strcpy(modelsav,model); 
Line 5328  int decodemodel ( char model[], int last Line 5396  int decodemodel ( char model[], int last
     return (1);      return (1);
 }  }
   
 calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )  int calandcheckages(int imx, int maxwav, double *agemin, double *agemax, int *nberr, int *nbwarn )
 {  {
   int i, m;    int i, m;
   
Line 5339  calandcheckages(int imx, int maxwav, dou Line 5407  calandcheckages(int imx, int maxwav, dou
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){        if((int)moisdc[i]==99 && (int)andc[i]==9999 && s[m][i]>nlstate){
         *nberr++;          *nberr = *nberr + 1;
         printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);          printf("Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
         fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased\n",(int)moisdc[i],(int)andc[i],num[i],i);          fprintf(ficlog,"Error! Date of death (month %2d and year %4d) of individual %ld on line %d was unknown, you must set an arbitrary year of death or he/she is skipped and results are biased (%d)\n",(int)moisdc[i],(int)andc[i],num[i],i, *nberr);
         s[m][i]=-1;          s[m][i]=-1;
       }        }
       if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){        if((int)moisdc[i]==99 && (int)andc[i]!=9999 && s[m][i]>nlstate){
         *nberr++;          (*nberr)++;
         printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]);           printf("Error! Month of death of individual %ld on line %d was unknown %2d, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,(int)moisdc[i]); 
         fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]);           fprintf(ficlog,"Error! Month of death of individual %ld on line %d was unknown %f, you should set it otherwise the information on the death is skipped and results are biased.\n",num[i],i,moisdc[i]); 
         s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */          s[m][i]=-1; /* We prefer to skip it (and to skip it in version 0.8a1 too */
Line 5358  calandcheckages(int imx, int maxwav, dou Line 5426  calandcheckages(int imx, int maxwav, dou
     for(m=firstpass; (m<= lastpass); m++){      for(m=firstpass; (m<= lastpass); m++){
       if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){        if(s[m][i] >0 || s[m][i]==-2 || s[m][i]==-4 || s[m][i]==-5){
         if (s[m][i] >= nlstate+1) {          if (s[m][i] >= nlstate+1) {
           if(agedc[i]>0)            if(agedc[i]>0){
             if((int)moisdc[i]!=99 && (int)andc[i]!=9999)              if((int)moisdc[i]!=99 && (int)andc[i]!=9999){
               agev[m][i]=agedc[i];                agev[m][i]=agedc[i];
           /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/            /*if(moisdc[i]==99 && andc[i]==9999) s[m][i]=-1;*/
             else {              }else {
               if ((int)andc[i]!=9999){                if ((int)andc[i]!=9999){
                 nbwarn++;                  nbwarn++;
                 printf("Warning negative age at death: %ld line:%d\n",num[i],i);                  printf("Warning negative age at death: %ld line:%d\n",num[i],i);
Line 5370  calandcheckages(int imx, int maxwav, dou Line 5438  calandcheckages(int imx, int maxwav, dou
                 agev[m][i]=-1;                  agev[m][i]=-1;
               }                }
             }              }
             } /* agedc > 0 */
         }          }
         else if(s[m][i] !=9){ /* Standard case, age in fractional          else if(s[m][i] !=9){ /* Standard case, age in fractional
                                  years but with the precision of a month */                                   years but with the precision of a month */
Line 5400  calandcheckages(int imx, int maxwav, dou Line 5469  calandcheckages(int imx, int maxwav, dou
   for (i=1; i<=imx; i++)  {    for (i=1; i<=imx; i++)  {
     for(m=firstpass; (m<=lastpass); m++){      for(m=firstpass; (m<=lastpass); m++){
       if (s[m][i] > (nlstate+ndeath)) {        if (s[m][i] > (nlstate+ndeath)) {
         *nberr++;          (*nberr)++;
         printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);               printf("Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
         fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);               fprintf(ficlog,"Error: on wave %d of individual %d status %d > (nlstate+ndeath)=(%d+%d)=%d\n",m,i,s[m][i],nlstate, ndeath, nlstate+ndeath);     
         return 1;          return 1;
Line 5425  calandcheckages(int imx, int maxwav, dou Line 5494  calandcheckages(int imx, int maxwav, dou
     return (1);      return (1);
 }  }
   
 syscompilerinfo()  #if defined(_MSC_VER)
   /*printf("Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/
   /*fprintf(ficlog, "Visual C++ compiler: %s \n;", _MSC_FULL_VER);*/
   //#include "stdafx.h"
   //#include <stdio.h>
   //#include <tchar.h>
   //#include <windows.h>
   //#include <iostream>
   typedef BOOL(WINAPI *LPFN_ISWOW64PROCESS) (HANDLE, PBOOL);
   
   LPFN_ISWOW64PROCESS fnIsWow64Process;
   
   BOOL IsWow64()
   {
           BOOL bIsWow64 = FALSE;
   
           //typedef BOOL (APIENTRY *LPFN_ISWOW64PROCESS)
           //  (HANDLE, PBOOL);
   
           //LPFN_ISWOW64PROCESS fnIsWow64Process;
   
           HMODULE module = GetModuleHandle(_T("kernel32"));
           const char funcName[] = "IsWow64Process";
           fnIsWow64Process = (LPFN_ISWOW64PROCESS)
                   GetProcAddress(module, funcName);
   
           if (NULL != fnIsWow64Process)
           {
                   if (!fnIsWow64Process(GetCurrentProcess(),
                           &bIsWow64))
                           //throw std::exception("Unknown error");
                           printf("Unknown error\n");
           }
           return bIsWow64 != FALSE;
   }
   #endif
   
   void syscompilerinfo()
  {   {
    /* #include "syscompilerinfo.h"*/     /* #include "syscompilerinfo.h"*/
 #include <gnu/libc-version.h>  
   #if defined __INTEL_COMPILER
 #if defined(__GNUC__)  #if defined(__GNUC__)
 # if defined(__GNUC_PATCHLEVEL__)          struct utsname sysInfo;  /* For Intel on Linux and OS/X */
 #  define __GNUC_VERSION__ (__GNUC__ * 10000 \  #endif
                             + __GNUC_MINOR__ * 100 \  #elif defined(__GNUC__) 
                             + __GNUC_PATCHLEVEL__)  #ifndef  __APPLE__
 # else  #include <gnu/libc-version.h>  /* Only on gnu */
 #  define __GNUC_VERSION__ (__GNUC__ * 10000 \  #endif
                             + __GNUC_MINOR__ * 100)     struct utsname sysInfo;
 # endif     int cross = CROSS;
      if (cross){
              printf("Cross-");
              fprintf(ficlog, "Cross-");
      }
 #endif  #endif
   
   #include <stdint.h>
   
      printf("Compiled with:");fprintf(ficlog,"Compiled with:");
   #if defined(__clang__)
      printf(" Clang/LLVM");fprintf(ficlog," Clang/LLVM"); /* Clang/LLVM. ---------------------------------------------- */
   #endif
   #if defined(__ICC) || defined(__INTEL_COMPILER)
      printf(" Intel ICC/ICPC");fprintf(ficlog," Intel ICC/ICPC");/* Intel ICC/ICPC. ------------------------------------------ */
   #endif
   #if defined(__GNUC__) || defined(__GNUG__)
      printf(" GNU GCC/G++");fprintf(ficlog," GNU GCC/G++");/* GNU GCC/G++. --------------------------------------------- */
   #endif
   #if defined(__HP_cc) || defined(__HP_aCC)
      printf(" Hewlett-Packard C/aC++");fprintf(fcilog," Hewlett-Packard C/aC++"); /* Hewlett-Packard C/aC++. ---------------------------------- */
   #endif
   #if defined(__IBMC__) || defined(__IBMCPP__)
      printf(" IBM XL C/C++"); fprintf(ficlog," IBM XL C/C++");/* IBM XL C/C++. -------------------------------------------- */
   #endif
   #if defined(_MSC_VER)
      printf(" Microsoft Visual Studio");fprintf(ficlog," Microsoft Visual Studio");/* Microsoft Visual Studio. --------------------------------- */
   #endif
   #if defined(__PGI)
      printf(" Portland Group PGCC/PGCPP");fprintf(ficlog," Portland Group PGCC/PGCPP");/* Portland Group PGCC/PGCPP. ------------------------------- */
   #endif
   #if defined(__SUNPRO_C) || defined(__SUNPRO_CC)
      printf(" Oracle Solaris Studio");fprintf(ficlog," Oracle Solaris Studio\n");/* Oracle Solaris Studio. ----------------------------------- */
   #endif
      printf(" for ");fprintf(ficlog," for ");
      
 // http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros  // http://stackoverflow.com/questions/4605842/how-to-identify-platform-compiler-from-preprocessor-macros
 #ifdef _WIN32 // note the underscore: without it, it's not msdn official!  #ifdef _WIN32 // note the underscore: without it, it's not msdn official!
     // Windows (x64 and x86)      // Windows (x64 and x86)
      printf("Windows (x64 and x86) ");fprintf(ficlog,"Windows (x64 and x86) ");
 #elif __unix__ // all unices, not all compilers  #elif __unix__ // all unices, not all compilers
     // Unix      // Unix
      printf("Unix ");fprintf(ficlog,"Unix ");
 #elif __linux__  #elif __linux__
     // linux      // linux
      printf("linux ");fprintf(ficlog,"linux ");
 #elif __APPLE__  #elif __APPLE__
     // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though...      // Mac OS, not sure if this is covered by __posix__ and/or __unix__ though..
      printf("Mac OS ");fprintf(ficlog,"Mac OS ");
 #endif  #endif
   
 /*  __MINGW32__   */  /*  __MINGW32__   */
Line 5460  syscompilerinfo() Line 5604  syscompilerinfo()
 /* _WIN64  // Defined for applications for Win64. */  /* _WIN64  // Defined for applications for Win64. */
 /* _M_X64 // Defined for compilations that target x64 processors. */  /* _M_X64 // Defined for compilations that target x64 processors. */
 /* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */  /* _DEBUG // Defined when you compile with /LDd, /MDd, and /MTd. */
 #include <stdint.h>  
 #if UINTPTR_MAX == 0xffffffff  #if UINTPTR_MAX == 0xffffffff
    printf("32-bit \n"); /* 32-bit */     printf(" 32-bit"); fprintf(ficlog," 32-bit");/* 32-bit */
 #elif UINTPTR_MAX == 0xffffffffffffffff  #elif UINTPTR_MAX == 0xffffffffffffffff
   printf("64-bit \n");/* 64-bit */     printf(" 64-bit"); fprintf(ficlog," 64-bit");/* 64-bit */
 #else  #else
  printf("wtf-bit \n"); /* wtf */     printf(" wtf-bit"); fprintf(ficlog," wtf-bit");/* wtf */
 #endif  #endif
   
 struct utsname sysInfo;  #if defined(__GNUC__)
   # if defined(__GNUC_PATCHLEVEL__)
   #  define __GNUC_VERSION__ (__GNUC__ * 10000 \
                               + __GNUC_MINOR__ * 100 \
                               + __GNUC_PATCHLEVEL__)
   # else
   #  define __GNUC_VERSION__ (__GNUC__ * 10000 \
                               + __GNUC_MINOR__ * 100)
   # endif
      printf(" using GNU C version %d.\n", __GNUC_VERSION__);
      fprintf(ficlog, " using GNU C version %d.\n", __GNUC_VERSION__);
   
    if (uname(&sysInfo) != -1) {     if (uname(&sysInfo) != -1) {
       puts(sysInfo.sysname);       printf("Running on: %s %s %s %s %s\n",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
       puts(sysInfo.nodename);       fprintf(ficlog,"Running on: %s %s %s %s %s\n ",sysInfo.sysname, sysInfo.nodename, sysInfo.release, sysInfo.version, sysInfo.machine);
       puts(sysInfo.release);  
       puts(sysInfo.version);  
       puts(sysInfo.machine);  
    }     }
    else     else
       perror("uname() error");        perror("uname() error");
    printf("GNU C version %d\n", __GNUC_VERSION__);     //#ifndef __INTEL_COMPILER 
   printf("GNU libc version: %s\n", gnu_get_libc_version());  #if !defined (__INTEL_COMPILER) && !defined(__APPLE__)
      printf("GNU libc version: %s\n", gnu_get_libc_version()); 
      fprintf(ficlog,"GNU libc version: %s\n", gnu_get_libc_version());
   #endif
   #endif
   
      //   void main()
      //   {
   #if defined(_MSC_VER)
      if (IsWow64()){
              printf("The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
              fprintf(ficlog, "The program (probably compiled for 32bit) is running under WOW64 (64bit) emulation.\n");
      }
      else{
              printf("The process is not running under WOW64 (i.e probably on a 64bit Windows).\n");
              fprintf(ficlog,"The programm is not running under WOW64 (i.e probably on a 64bit Windows).\n");
      }
      //      printf("\nPress Enter to continue...");
      //      getchar();
      //   }
   
   #endif
      
   
  }   }
   
   int prevalence_limit(double *p, double **prlim, double ageminpar, double agemaxpar){
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
     int i, j, k, i1 ;
     double ftolpl = 1.e-10;
     double age, agebase, agelim;
   
       strcpy(filerespl,"pl");
       strcat(filerespl,fileres);
       if((ficrespl=fopen(filerespl,"w"))==NULL) {
         printf("Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
         fprintf(ficlog,"Problem with period (stable) prevalence resultfile: %s\n", filerespl);return 1;
       }
       printf("Computing period (stable) prevalence: result on file '%s' \n", filerespl);
       fprintf(ficlog,"Computing period (stable) prevalence: result on file '%s' \n", filerespl);
       pstamp(ficrespl);
       fprintf(ficrespl,"# Period (stable) prevalence \n");
       fprintf(ficrespl,"#Age ");
       for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
       fprintf(ficrespl,"\n");
     
       /* prlim=matrix(1,nlstate,1,nlstate);*/ /* back in main */
   
       agebase=ageminpar;
       agelim=agemaxpar;
   
       i1=pow(2,cptcoveff);
       if (cptcovn < 1){i1=1;}
   
       for(cptcov=1,k=0;cptcov<=i1;cptcov++){
       /* for(cptcov=1,k=0;cptcov<=1;cptcov++){ */
         //for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
           k=k+1;
           /* to clean */
           //printf("cptcov=%d cptcod=%d codtab=%d\n",cptcov, cptcod,codtab[cptcod][cptcov]);
           fprintf(ficrespl,"\n#******");
           printf("\n#******");
           fprintf(ficlog,"\n#******");
           for(j=1;j<=cptcoveff;j++) {
             fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             printf(" V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             fprintf(ficlog," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           }
           fprintf(ficrespl,"******\n");
           printf("******\n");
           fprintf(ficlog,"******\n");
   
           fprintf(ficrespl,"#Age ");
           for(j=1;j<=cptcoveff;j++) {
             fprintf(ficrespl,"V%d %d",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           }
           for(i=1; i<=nlstate;i++) fprintf(ficrespl,"%d-%d ",i,i);
           fprintf(ficrespl,"\n");
           
           for (age=agebase; age<=agelim; age++){
           /* for (age=agebase; age<=agebase; age++){ */
             prevalim(prlim, nlstate, p, age, oldm, savm,ftolpl,k);
             fprintf(ficrespl,"%.0f ",age );
             for(j=1;j<=cptcoveff;j++)
               fprintf(ficrespl,"%d %d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
             for(i=1; i<=nlstate;i++)
               fprintf(ficrespl," %.5f", prlim[i][i]);
             fprintf(ficrespl,"\n");
           } /* Age */
           /* was end of cptcod */
       } /* cptcov */
   }
   
   int hPijx(double *p, int bage, int fage){
       /*------------- h Pij x at various ages ------------*/
   
     int stepsize;
     int agelim;
     int hstepm;
     int nhstepm;
     int h, i, i1, j, k;
   
     double agedeb;
     double ***p3mat;
   
       strcpy(filerespij,"pij");  strcat(filerespij,fileres);
       if((ficrespij=fopen(filerespij,"w"))==NULL) {
         printf("Problem with Pij resultfile: %s\n", filerespij); return 1;
         fprintf(ficlog,"Problem with Pij resultfile: %s\n", filerespij); return 1;
       }
       printf("Computing pij: result on file '%s' \n", filerespij);
       fprintf(ficlog,"Computing pij: result on file '%s' \n", filerespij);
     
       stepsize=(int) (stepm+YEARM-1)/YEARM;
       /*if (stepm<=24) stepsize=2;*/
   
       agelim=AGESUP;
       hstepm=stepsize*YEARM; /* Every year of age */
       hstepm=hstepm/stepm; /* Typically 2 years, = 2/6 months = 4 */ 
   
       /* hstepm=1;   aff par mois*/
       pstamp(ficrespij);
       fprintf(ficrespij,"#****** h Pij x Probability to be in state j at age x+h being in i at x ");
       i1= pow(2,cptcoveff);
      for(cptcov=1,k=0;cptcov<=i1;cptcov++){
         /*for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
           k=k+1; 
       /* for (k=1; k <= (int) pow(2,cptcoveff); k++){*/
           fprintf(ficrespij,"\n#****** ");
           for(j=1;j<=cptcoveff;j++) 
             fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
           fprintf(ficrespij,"******\n");
           
           for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
             nhstepm=(int) rint((agelim-agedeb)*YEARM/stepm); /* Typically 20 years = 20*12/6=40 */ 
             nhstepm = nhstepm/hstepm; /* Typically 40/4=10 */
   
             /*      nhstepm=nhstepm*YEARM; aff par mois*/
   
             p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             oldm=oldms;savm=savms;
             hpxij(p3mat,nhstepm,agedeb,hstepm,p,nlstate,stepm,oldm,savm, k);  
             fprintf(ficrespij,"# Cov Agex agex+h hpijx with i,j=");
             for(i=1; i<=nlstate;i++)
               for(j=1; j<=nlstate+ndeath;j++)
                 fprintf(ficrespij," %1d-%1d",i,j);
             fprintf(ficrespij,"\n");
             for (h=0; h<=nhstepm; h++){
               /*agedebphstep = agedeb + h*hstepm/YEARM*stepm;*/
               fprintf(ficrespij,"%d %3.f %3.f",k, agedeb, agedeb + h*hstepm/YEARM*stepm );
               for(i=1; i<=nlstate;i++)
                 for(j=1; j<=nlstate+ndeath;j++)
                   fprintf(ficrespij," %.5f", p3mat[i][j][h]);
               fprintf(ficrespij,"\n");
             }
             free_ma3x(p3mat,1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
             fprintf(ficrespij,"\n");
           }
         /*}*/
       }
   }
   
   
 /***********************************************/  /***********************************************/
 /**************** Main Program *****************/  /**************** Main Program *****************/
 /***********************************************/  /***********************************************/
Line 5634  int main(int argc, char *argv[]) Line 5944  int main(int argc, char *argv[])
   strcpy(command,"mkdir ");    strcpy(command,"mkdir ");
   strcat(command,optionfilefiname);    strcat(command,optionfilefiname);
   if((outcmd=system(command)) != 0){    if((outcmd=system(command)) != 0){
     printf("Problem creating directory or it already exists %s%s, err=%d\n",path,optionfilefiname,outcmd);      printf("Directory already exists (or can't create it) %s%s, err=%d\n",path,optionfilefiname,outcmd);
     /* fprintf(ficlog,"Problem creating directory %s%s\n",path,optionfilefiname); */      /* fprintf(ficlog,"Problem creating directory %s%s\n",path,optionfilefiname); */
     /* fclose(ficlog); */      /* fclose(ficlog); */
 /*     exit(1); */  /*     exit(1); */
Line 6614  Interval (in months) between two waves: Line 6924  Interval (in months) between two waves:
   
   
     /*--------------- Prevalence limit  (period or stable prevalence) --------------*/      /*--------------- Prevalence limit  (period or stable prevalence) --------------*/
 #include "prevlim.h"  /* Use ficrespl, ficlog */      /*#include "prevlim.h"*/  /* Use ficrespl, ficlog */
       prlim=matrix(1,nlstate,1,nlstate);
       prevalence_limit(p, prlim,  ageminpar, agemaxpar);
     fclose(ficrespl);      fclose(ficrespl);
   
 #ifdef FREEEXIT2  #ifdef FREEEXIT2
Line 6622  Interval (in months) between two waves: Line 6934  Interval (in months) between two waves:
 #endif  #endif
   
     /*------------- h Pij x at various ages ------------*/      /*------------- h Pij x at various ages ------------*/
 #include "hpijx.h"      /*#include "hpijx.h"*/
       hPijx(p, bage, fage);
     fclose(ficrespij);      fclose(ficrespij);
   
   /*-------------- Variance of one-step probabilities---*/    /*-------------- Variance of one-step probabilities---*/
Removed from v.1.167  
changed lines
  Added in v.1.181

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