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

Diff for /imach/src/imach.c between versions 1.268 and 1.275

version 1.268, 2017/05/18 20:09:32 version 1.275, 2017/06/30 13:39:33
Line 1 Line 1
 /* $Id$  /* $Id$
   $State$    $State$
   $Log$    $Log$
     Revision 1.275  2017/06/30 13:39:33  brouard
     Summary: Saito's color
   
     Revision 1.274  2017/06/29 09:47:08  brouard
     Summary: Version 0.99r14
   
     Revision 1.273  2017/06/27 11:06:02  brouard
     Summary: More documentation on projections
   
     Revision 1.272  2017/06/27 10:22:40  brouard
     Summary: Color of backprojection changed from 6 to 5(yellow)
   
     Revision 1.271  2017/06/27 10:17:50  brouard
     Summary: Some bug with rint
   
     Revision 1.270  2017/05/24 05:45:29  brouard
     *** empty log message ***
   
     Revision 1.269  2017/05/23 08:39:25  brouard
     Summary: Code into subroutine, cleanings
   
   Revision 1.268  2017/05/18 20:09:32  brouard    Revision 1.268  2017/05/18 20:09:32  brouard
   Summary: backprojection and confidence intervals of backprevalence    Summary: backprojection and confidence intervals of backprevalence
   
Line 1084  FILE *ficrescveij; Line 1105  FILE *ficrescveij;
 char filerescve[FILENAMELENGTH];  char filerescve[FILENAMELENGTH];
 FILE  *ficresvij;  FILE  *ficresvij;
 char fileresv[FILENAMELENGTH];  char fileresv[FILENAMELENGTH];
 FILE  *ficresvpl;  
 char fileresvpl[FILENAMELENGTH];  
 FILE  *ficresvbl;  
 char fileresvbl[FILENAMELENGTH];  
 char title[MAXLINE];  char title[MAXLINE];
 char model[MAXLINE]; /**< The model line */  char model[MAXLINE]; /**< The model line */
 char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH],  fileresplb[FILENAMELENGTH];  char optionfile[FILENAMELENGTH], datafile[FILENAMELENGTH],  filerespl[FILENAMELENGTH],  fileresplb[FILENAMELENGTH];
Line 2932  double **pmij(double **ps, double *cov, Line 2950  double **pmij(double **ps, double *cov,
   /* Diag(w_x) */    /* Diag(w_x) */
   /* Problem with prevacurrent which can be zero */    /* Problem with prevacurrent which can be zero */
   sumnew=0.;    sumnew=0.;
   /* for (ii=1;ii<=nlstate+ndeath;ii++){ */    /*for (ii=1;ii<=nlstate+ndeath;ii++){*/
   for (ii=1;ii<=nlstate;ii++){ /* Only on live states */    for (ii=1;ii<=nlstate;ii++){ /* Only on live states */
       /* printf(" agefin=%d, ii=%d, ij=%d, prev=%f\n",(int)agefin,ii, ij, prevacurrent[(int)agefin][ii][ij]);  */
     sumnew+=prevacurrent[(int)agefin][ii][ij];      sumnew+=prevacurrent[(int)agefin][ii][ij];
   }    }
   if(sumnew >0.01){  /* At least some value in the prevalence */    if(sumnew >0.01){  /* At least some value in the prevalence */
     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++)
       doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);          doldm[ii][j]=(ii==j ? prevacurrent[(int)agefin][ii][ij]/sumnew : 0.0);
     }      }
   }else{    }else{
     for (ii=1;ii<=nlstate+ndeath;ii++){      for (ii=1;ii<=nlstate+ndeath;ii++){
Line 3239  double ***hbxij(double ***po, int nhstep Line 3258  double ***hbxij(double ***po, int nhstep
       newm=savm;        newm=savm;
       /* Covariates have to be included here again */        /* Covariates have to be included here again */
       cov[1]=1.;        cov[1]=1.;
       agexact=age-((h-1)*hstepm + (d))*stepm/YEARM; /* age just before transition, d or d-1? */        agexact=age-( (h-1)*hstepm + (d)  )*stepm/YEARM; /* age just before transition, d or d-1? */
       /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */        /* agexact=age+((h-1)*hstepm + (d-1))*stepm/YEARM; /\* age just before transition *\/ */
       cov[2]=agexact;        cov[2]=agexact;
       if(nagesqr==1)        if(nagesqr==1)
Line 3841  void likelione(FILE *ficres,double p[], Line 3860  void likelione(FILE *ficres,double p[],
     else if(mle >=1)      else if(mle >=1)
       fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);        fprintf(fichtm,"\n<br>File of contributions to the likelihood computed with optimized parameters mle = %d.",mle);
     fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));      fprintf(fichtm," You should at least run with mle >= 1 to get starting values corresponding to the optimized parameters in order to visualize the real contribution of each individual/wave: <a href=\"%s\">%s</a><br>\n",subdirf(fileresilk),subdirf(fileresilk));
           fprintf(fichtm,"\n<br>Equation of the model: <b>model=1+age+%s</b><br>\n",model); 
               
     for (k=1; k<= nlstate ; k++) {      for (k=1; k<= nlstate ; k++) {
       fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \        fprintf(fichtm,"<br>- Probability p<sub>%dj</sub> by origin %d and destination j. Dot's sizes are related to corresponding weight: <a href=\"%s-p%dj.png\">%s-p%dj.png</a><br> \
Line 5463  void  concatwav(int wav[], int **dh, int Line 5482  void  concatwav(int wav[], int **dh, int
   /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm.     /* hstepm beeing the number of stepms, if hstepm=1 the length of hstepm is stepm. 
      nhstepm is the number of hstepm from age to agelim        nhstepm is the number of hstepm from age to agelim 
      nstepm is the number of stepm from age to agelin.        nstepm is the number of stepm from age to agelin. 
      Look at hpijx to understand the reason of that which relies in memory size       Look at hpijx to understand the reason which relies in memory size consideration
      and note for a fixed period like estepm months */       and note for a fixed period like estepm months */
   /* We decided (b) to get a life expectancy respecting the most precise curvature of the    /* We decided (b) to get a life expectancy respecting the most precise curvature of the
      survival function given by stepm (the optimization length). Unfortunately it       survival function given by stepm (the optimization length). Unfortunately it
Line 5694  void  concatwav(int wav[], int **dh, int Line 5713  void  concatwav(int wav[], int **dh, int
           /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/            /* if((int)age==70)printf("i=%2d,j=%2d,h=%2d,age=%3d,%9.4f,%9.4f,%9.4f\n",i,j,h,(int)age,p3mat[i][j][h],hf,eij[i][j][(int)age]);*/
                                                                                   
         }          }
                   
       /* Standard deviation of expectancies ij */         
     fprintf(ficresstdeij,"%3.0f",age );      fprintf(ficresstdeij,"%3.0f",age );
     for(i=1; i<=nlstate;i++){      for(i=1; i<=nlstate;i++){
       eip=0.;        eip=0.;
Line 5709  void  concatwav(int wav[], int **dh, int Line 5729  void  concatwav(int wav[], int **dh, int
     }      }
     fprintf(ficresstdeij,"\n");      fprintf(ficresstdeij,"\n");
                                   
       /* Variance of expectancies ij */           
     fprintf(ficrescveij,"%3.0f",age );      fprintf(ficrescveij,"%3.0f",age );
     for(i=1; i<=nlstate;i++)      for(i=1; i<=nlstate;i++)
       for(j=1; j<=nlstate;j++){        for(j=1; j<=nlstate;j++){
Line 6062  void  concatwav(int wav[], int **dh, int Line 6083  void  concatwav(int wav[], int **dh, int
  }  /* end varevsij */   }  /* end varevsij */
   
 /************ Variance of prevlim ******************/  /************ Variance of prevlim ******************/
  void varprevlim(char fileres[], double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres)   void varprevlim(char fileresvpl[], FILE *ficresvpl, double **varpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **prlim, double ftolpl, int *ncvyearp, int ij, char strstart[], int nres)
 {  {
   /* Variance of prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/    /* Variance of prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
   /*  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);*/
Line 6187  void  concatwav(int wav[], int **dh, int Line 6208  void  concatwav(int wav[], int **dh, int
   
   
 /************ Variance of backprevalence limit ******************/  /************ Variance of backprevalence limit ******************/
  void varbrevlim(char fileres[], double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres)   void varbrevlim(char fileresvbl[], FILE  *ficresvbl, double **varbpl, double **matcov, double x[], double delti[], int nlstate, int stepm, double bage, double fage, double **oldm, double **savm, double **bprlim, double ftolpl, int mobilavproj, int *ncvyearp, int ij, char strstart[], int nres)
 {  {
   /* Variance of backward prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/    /* Variance of backward prevalence limit  for each state ij using current parameters x[] and estimates of neighbourhood give by delti*/
   /*  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);*/
Line 6652  void printinghtml(char fileresu[], char Line 6673  void printinghtml(char fileresu[], char
                   int lastpass, int stepm, int weightopt, char model[],\                    int lastpass, int stepm, int weightopt, char model[],\
                   int imx,int jmin, int jmax, double jmeanint,char rfileres[],\                    int imx,int jmin, int jmax, double jmeanint,char rfileres[],\
                   int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \                    int popforecast, int mobilav, int prevfcast, int mobilavproj, int backcast, int estepm , \
                   double jprev1, double mprev1,double anprev1, double dateprev1, \                    double jprev1, double mprev1,double anprev1, double dateprev1, double dateproj1, double dateback1, \
                   double jprev2, double mprev2,double anprev2, double dateprev2){                    double jprev2, double mprev2,double anprev2, double dateprev2, double dateproj2, double dateback2){
   int jj1, k1, i1, cpt, k4, nres;    int jj1, k1, i1, cpt, k4, nres;
   
    fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \     fprintf(fichtm,"<ul><li><a href='#firstorder'>Result files (first order: no variance)</a>\n \
Line 6807  divided by h: <sub>h</sub>P<sub>ij</sub> Line 6828  divided by h: <sub>h</sub>P<sub>ij</sub>
      if(prevfcast==1){       if(prevfcast==1){
        /* Projection of prevalence up to period (stable) prevalence in each health state */         /* Projection of prevalence up to period (stable) prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){         for(cpt=1; cpt<=nlstate;cpt++){
          fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \           fprintf(fichtm,"<br>\n- Projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), from year %.1f up to year %.1f tending to period (stable) prevalence in state %d. Or probability to be in state %d being in an observed weighted state (from 1 to %d). <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
 <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);  <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateproj1, dateproj2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJ_"),cpt,k1,nres);
        }         }
      }       }
      if(backcast==1){       if(backcast==1){
       /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */        /* Back projection of prevalence up to stable (mixed) back-prevalence in each health state */
        for(cpt=1; cpt<=nlstate;cpt++){         for(cpt=1; cpt<=nlstate;cpt++){
          fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d) up to stable (mixed) back prevalence in state %d. Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) with weights corresponding to observed prevalence  at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \           fprintf(fichtm,"<br>\n- Back projection of cross-sectional prevalence (estimated with cases observed from %.1f to %.1f and mobil_average=%d), \
 <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);   from year %.1f up to year %.1f (probably close to stable [mixed] back prevalence in state %d (randomness in cross-sectional prevalence is not taken into \
    account but can visually be appreciated). Or probability to have been in an state %d, knowing that the person was in either state (1 or %d) \
   with weights corresponding to observed prevalence at different ages. <a href=\"%s_%d-%d-%d.svg\">%s_%d-%d-%d.svg</a><br> \
    <img src=\"%s_%d-%d-%d.svg\">", dateprev1, dateprev2, mobilavproj, dateback1, dateback2, cpt, cpt, nlstate, subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres,subdirf2(optionfilefiname,"PROJB_"),cpt,k1,nres);
        }         }
      }       }
                     
Line 6921  true period expectancies (those weighted Line 6945  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[], int offyear, int offbyear){  void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int backcast, char pathc[], double p[], int offyear, int offbyear){
   
   char dirfileres[132],optfileres[132];    char dirfileres[132],optfileres[132];
   char gplotcondition[132], gplotlabel[132];    char gplotcondition[132], gplotlabel[132];
Line 6930  void printinggnuplot(char fileresu[], ch Line 6954  void printinggnuplot(char fileresu[], ch
   int ng=0;    int ng=0;
   int vpopbased;    int vpopbased;
   int ioffset; /* variable offset for columns */    int ioffset; /* variable offset for columns */
     int iyearc=1; /* variable column for year of projection  */
     int iagec=1; /* variable column for age of projection  */
   int nres=0; /* Index of resultline */    int nres=0; /* Index of resultline */
   int istart=1; /* For starting graphs in projections */    int istart=1; /* For starting graphs in projections */
   
Line 6943  void printinggnuplot(char fileresu[], ch Line 6969  void printinggnuplot(char fileresu[], ch
   /*#endif */    /*#endif */
   m=pow(2,cptcoveff);    m=pow(2,cptcoveff);
   
     /* diagram of the model */
     fprintf(ficgp,"\n#Diagram of the model \n");
     fprintf(ficgp,"\ndelta=0.03;delta2=0.07;unset arrow;\n");
     fprintf(ficgp,"yoff=(%d > 2? 0:1);\n",nlstate);
     fprintf(ficgp,"\n#Peripheral arrows\nset for [i=1:%d] for [j=1:%d] arrow i*10+j from cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.95*(cos(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0) - cos(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta2:0)), -0.95*(sin(pi*((1-(%d/2)*2./%d)/2+(i-1)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) - sin(pi*((1-(%d/2)*2./%d)/2+(j-1)*2./%d))+( i!=j?(i-j)/abs(i-j)*delta2:0)) ls (i < j? 1:2)\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
   
     fprintf(ficgp,"\n#Centripete arrows (turning in other direction (1-i) instead of (i-1)) \nset for [i=1:%d] arrow (%d+1)*10+i from cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))-(i!=j?(i-j)/abs(i-j)*delta:0), yoff +sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) rto -0.80*(cos(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d))+(i!=j?(i-j)/abs(i-j)*delta:0)  ), -0.80*(sin(pi*((1-(%d/2)*2./%d)/2+(1-i)*2./%d)) + (i!=j?(i-j)/abs(i-j)*delta:0) + yoff ) ls 4\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
     fprintf(ficgp,"\n#show arrow\nunset label\n");
     fprintf(ficgp,"\n#States labels, starting from 2 (2-i) instead of (1-i), was (i-1)\nset for [i=1:%d] label i sprintf(\"State %%d\",i) center at cos(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)), yoff+sin(pi*((1-(%d/2)*2./%d)/2+(2-i)*2./%d)) font \"helvetica, 16\" tc rgbcolor \"blue\"\n",nlstate,nlstate,nlstate,nlstate,nlstate,nlstate,nlstate);
     fprintf(ficgp,"\nset label %d+1 sprintf(\"State %%d\",%d+1) center at 0.,0.  font \"helvetica, 16\" tc rgbcolor \"red\"\n",nlstate,nlstate);
     fprintf(ficgp,"\n#show label\nunset border;unset xtics; unset ytics;\n");
     fprintf(ficgp,"\n\nset ter svg size 640, 480;set out \"%s_.svg\" \n",subdirf2(optionfilefiname,"D_"));
     fprintf(ficgp,"unset log y; plot [-1.2:1.2][yoff-1.2:1.2] 1/0 not; set out;reset;\n");
   
   /* Contribution to likelihood */    /* Contribution to likelihood */
   /* Plot the probability implied in the likelihood */    /* Plot the probability implied in the likelihood */
   fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");    fprintf(ficgp,"\n# Contributions to the Likelihood, mle >=1. For mle=4 no interpolation, pure matrix products.\n#\n");
Line 7034  void printinggnuplot(char fileresu[], ch Line 7074  void printinggnuplot(char fileresu[], ch
                   
         fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_"));          fprintf(ficgp,"\" t\"\" w l lt 1,\"%s\" u 1:((",subdirf2(fileresu,"P_"));
         if(cptcoveff ==0){          if(cptcoveff ==0){
           fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3",      2+(cpt-1),  cpt );            fprintf(ficgp,"$%d)) t 'Observed prevalence in state %d' with line lt 3",      2+3*(cpt-1),  cpt );
         }else{          }else{
           kl=0;            kl=0;
           for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */            for (k=1; k<=cptcoveff; k++){    /* For each combination of covariate  */
Line 7092  void printinggnuplot(char fileresu[], ch Line 7132  void printinggnuplot(char fileresu[], ch
               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\"Backward (stable) prevalence\" w l lt 0,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);              fprintf(ficgp,"\" t\"Backward (stable) prevalence\" w l lt 6 dt 3,\"%s\" every :::%d::%d u 1:($2==%d ? $3+1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);
             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==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres);               fprintf(ficgp,"\" t\"95%% CI\" w l lt 5,\"%s\" every :::%d::%d u 1:($2==%d ? $3-1.96*$4 : 1/0) \"%%lf %%lf",subdirf2(fileresu,"VBL_"),nres-1,nres-1,nres); 
             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");              fprintf(ficgp,"\" t\"\" w l lt 4");
           } /* end if backprojcast */            } /* end if backprojcast */
         } /* end if backcast */          } /* end if backcast */
         fprintf(ficgp,"\nset out ;unset label;\n");          fprintf(ficgp,"\nset out ;unset label;\n");
Line 7503  set ter svg size 640, 480\nunset log y\n Line 7543  set ter svg size 640, 480\nunset log y\n
             /*#  1    2        3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */              /*#  1    2        3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
             fprintf(ficgp," u %d:(", ioffset);               fprintf(ficgp," u %d:(", ioffset); 
             if(i==nlstate+1){              if(i==nlstate+1){
               fprintf(ficgp," $%d/(1.-$%d)):5 t 'pw.%d' with line lc variable ",        \                fprintf(ficgp," $%d/(1.-$%d)):1 t 'pw.%d' with line lc variable ",        \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
               fprintf(ficgp,",\\\n '' ");                fprintf(ficgp,",\\\n '' ");
               fprintf(ficgp," u %d:(",ioffset);                 fprintf(ficgp," u %d:(",ioffset); 
               fprintf(ficgp," (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", \                fprintf(ficgp," (($1-$2) == %d ) ? $%d/(1.-$%d) : 1/0):1 with labels center not ", \
                      offyear,                           \                       offyear,                           \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate );
             }else              }else
               fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ",      \                fprintf(ficgp," $%d/(1.-$%d)) t 'p%d%d' with line ",      \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,i,cpt );
           }else{ /* more than 2 covariates */            }else{ /* more than 2 covariates */
             if(cptcoveff ==1){              ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/
               ioffset=4; /* Age is in 4 */              /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
             }else{              /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
               ioffset=6; /* Age is in 6 */              iyearc=ioffset-1;
               /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/              iagec=ioffset;
               /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */  
             }     
             fprintf(ficgp," u %d:(",ioffset);               fprintf(ficgp," u %d:(",ioffset); 
             kl=0;              kl=0;
             strcpy(gplotcondition,"(");              strcpy(gplotcondition,"(");
Line 7542  set ter svg size 640, 480\nunset log y\n Line 7580  set ter svg size 640, 480\nunset log y\n
             /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */               /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ 
             /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/              /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
             if(i==nlstate+1){              if(i==nlstate+1){
               fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):5 t 'p.%d' with line lc variable", gplotcondition, \                fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0):%d t 'p.%d' with line lc variable", gplotcondition, \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,iyearc, cpt );
               fprintf(ficgp,",\\\n '' ");                fprintf(ficgp,",\\\n '' ");
               fprintf(ficgp," u %d:(",ioffset);                 fprintf(ficgp," u %d:(",iagec); 
               fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \                fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d/(1.-$%d) : 1/0):%d with labels center not ", gplotcondition, \
                      offyear,                           \                        iyearc, iagec, offyear,                           \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate, iyearc );
 /*  '' u 6:(($1==1 && $2==0  && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/  /*  '' u 6:(($1==1 && $2==0  && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
             }else{              }else{
               fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \                fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \
Line 7618  set ter svg size 640, 480\nunset log y\n Line 7656  set ter svg size 640, 480\nunset log y\n
             /*#  1    2        3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */              /*#  1    2        3   4   5  6    7  8   9   10  11  12  13  14  15  16  17  18 */
             fprintf(ficgp," u %d:(", ioffset);               fprintf(ficgp," u %d:(", ioffset); 
             if(i==nlstate+1){              if(i==nlstate+1){
               fprintf(ficgp," $%d/(1.-$%d)):5 t 'bw%d' with line lc variable ", \                fprintf(ficgp," $%d/(1.-$%d)):1 t 'bw%d' with line lc variable ", \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt );
               fprintf(ficgp,",\\\n '' ");                fprintf(ficgp,",\\\n '' ");
               fprintf(ficgp," u %d:(",ioffset);                 fprintf(ficgp," u %d:(",ioffset); 
               fprintf(ficgp," (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", \                fprintf(ficgp," (($1-$2) == %d ) ? $%d : 1/0):1 with labels center not ", \
                      offbyear,                          \                       offbyear,                          \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );
             }else              }else
               fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ",      \                fprintf(ficgp," $%d/(1.-$%d)) t 'b%d%d' with line ",      \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),  ioffset+1+(i-1)+(nlstate+1)*nlstate,cpt,i );
           }else{ /* more than 2 covariates */            }else{ /* more than 2 covariates */
             if(cptcoveff ==1){              ioffset=2*cptcoveff+2; /* Age is in 4 or 6 or etc.*/
               ioffset=4; /* Age is in 4 */              /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/
             }else{              /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */
               ioffset=6; /* Age is in 6 */              iyearc=ioffset-1;
               /*#  V1  = 1  V2 =  0 yearproj age p11 p21 p.1 p12 p22 p.2 p13 p23 p.3*/              iagec=ioffset;
               /*#   1    2   3    4    5      6  7   8   9   10   11 12  13   14  15 */  
             }     
             fprintf(ficgp," u %d:(",ioffset);               fprintf(ficgp," u %d:(",ioffset); 
             kl=0;              kl=0;
             strcpy(gplotcondition,"(");              strcpy(gplotcondition,"(");
Line 7657  set ter svg size 640, 480\nunset log y\n Line 7693  set ter svg size 640, 480\nunset log y\n
             /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */               /*6+1+(i-1)+(nlstate+1)*nlstate; 6+1+(1-1) +(2+1)*2=13 */ 
             /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/              /* ''  u 6:(($1==1 && $2==0 && $3==2 && $4==0)? $9/(1.-$15) : 1/0):($5==2000? 3:2) t 'p.1' with line lc variable*/
             if(i==nlstate+1){              if(i==nlstate+1){
               fprintf(ficgp,"%s ? $%d : 1/0):5 t 'bw%d' with line lc variable", gplotcondition, \                fprintf(ficgp,"%s ? $%d : 1/0):%d t 'bw%d' with line lc variable", gplotcondition, \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1),cpt );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1),iyearc,cpt );
               fprintf(ficgp,",\\\n '' ");                fprintf(ficgp,",\\\n '' ");
               fprintf(ficgp," u %d:(",ioffset);                 fprintf(ficgp," u %d:(",iagec); 
               /* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */                /* fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d/(1.-$%d) : 1/0):5 with labels center not ", gplotcondition, \ */
               fprintf(ficgp,"%s && (($5-$6) == %d ) ? $%d : 1/0):5 with labels center not ", gplotcondition, \                fprintf(ficgp,"%s && (($%d-$%d) == %d ) ? $%d : 1/0):%d with labels center not ", gplotcondition, \
                      offbyear,                          \                        iyearc,iagec,offbyear,                            \
                       ioffset+(cpt-1)*(nlstate+1)+1+(i-1) );                        ioffset+(cpt-1)*(nlstate+1)+1+(i-1), iyearc );
 /*  '' u 6:(($1==1 && $2==0  && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/  /*  '' u 6:(($1==1 && $2==0  && $3==2 && $4==0) && (($5-$6) == 1947) ? $10/(1.-$22) : 1/0):5 with labels center boxed not*/
             }else{              }else{
               /* fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ */                /* fprintf(ficgp,"%s ? $%d/(1.-$%d) : 1/0) t 'p%d%d' with line ", gplotcondition, \ */
Line 7860  set ter svg size 640, 480\nunset log y\n Line 7896  set ter svg size 640, 480\nunset log y\n
             }              }
             fprintf(ficgp,")");              fprintf(ficgp,")");
             if(ng ==2)              if(ng ==2)
               fprintf(ficgp," t \"p%d%d\" ", k2,k);                fprintf(ficgp," w l lt (%d*%d+%d)%%%d+1 dt %d t \"p%d%d\" ", nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
             else /* ng= 3 */              else /* ng= 3 */
               fprintf(ficgp," t \"i%d%d\" ", k2,k);                fprintf(ficgp," w l lt (%d*%d+%d)%%%d+1 dt %d t \"i%d%d\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
           }else{ /* end ng <> 1 */            }else{ /* end ng <> 1 */
             if( k !=k2) /* logit p11 is hard to draw */              if( k !=k2) /* logit p11 is hard to draw */
               fprintf(ficgp," t \"logit(p%d%d)\" ", k2,k);                fprintf(ficgp," w l lt (%d*%d+%d)%%%d+1 dt %d t \"logit(p%d%d)\" ",  nlstate+ndeath, k2, k, nlstate+ndeath, k2, k2,k);
           }            }
           if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)            if ((k+k2)!= (nlstate*2+ndeath) && ng != 1)
             fprintf(ficgp,",");              fprintf(ficgp,",");
Line 8112  set ter svg size 640, 480\nunset log y\n Line 8148  set ter svg size 640, 480\nunset log y\n
     
   
 /************** Forecasting ******************/  /************** Forecasting ******************/
 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){   void prevforecast(char fileres[], double anproj1, double mproj1, double jproj1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double ***prev, 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 8151  void prevforecast(char fileres[], double Line 8187  void prevforecast(char fileres[], double
   if(estepm < stepm){    if(estepm < stepm){
     printf ("Problem %d lower than %d\n",estepm, stepm);      printf ("Problem %d lower than %d\n",estepm, stepm);
   }    }
   else  hstepm=estepm;       else{
       hstepm=estepm;   
     }
     if(estepm > stepm){ /* Yes every two year */
       stepsize=2;
     }
   
   hstepm=hstepm/stepm;     hstepm=hstepm/stepm; 
   yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and    yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
Line 8196  void prevforecast(char fileres[], double Line 8237  void prevforecast(char fileres[], double
     for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {      for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {
       fprintf(ficresf,"\n");        fprintf(ficresf,"\n");
       fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);           fprintf(ficresf,"\n# Forecasting at date %.lf/%.lf/%.lf ",jproj1,mproj1,anproj1+yearp);   
       for (agec=fage; agec>=(ageminpar-1); agec--){         /* for (agec=fage; agec>=(ageminpar-1); agec--){  */
         for (agec=fage; agec>=(bage); agec--){ 
         nhstepm=(int) rint((agelim-agec)*YEARM/stepm);           nhstepm=(int) rint((agelim-agec)*YEARM/stepm); 
         nhstepm = nhstepm/hstepm;           nhstepm = nhstepm/hstepm; 
         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
Line 8218  void prevforecast(char fileres[], double Line 8260  void prevforecast(char fileres[], double
           ppij=0.;            ppij=0.;
           for(i=1; i<=nlstate;i++) {            for(i=1; i<=nlstate;i++) {
             /* if (mobilav>=1)  */              /* if (mobilav>=1)  */
             ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][i][k];              ppij=ppij+p3mat[i][j][h]*prev[(int)agec][i][k];
             /* else { */ /* even if mobilav==-1 we use mobaverage */              /* else { */ /* even if mobilav==-1 we use mobaverage */
             /*  ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */              /*  ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
             /* } */              /* } */
Line 8239  void prevforecast(char fileres[], double Line 8281  void prevforecast(char fileres[], double
   
 }  }
   
 /* /\************** Back Forecasting ******************\/ */  /************** Back Forecasting ******************/
 void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){   void prevbackforecast(char fileres[], double ***prevacurrent, double anback1, double mback1, double jback1, double ageminpar, double agemax, double dateprev1, double dateprev2, int mobilav, double bage, double fage, int firstpass, int lastpass, double anback2, double p[], int cptcoveff){
   /* back1, year, month, day of starting backection    /* back1, year, month, day of starting backection
      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
      anback2 year of en of backection (same day and month as back1).       anback2 year of end of backprojection (same day and month as back1).
        prevacurrent and prev are prevalences.
   */    */
   int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;    int yearp, stepsize, hstepm, nhstepm, j, k, cptcod, i, h, i1, k4, nres=0;
   double agec; /* generic age */    double agec; /* generic age */
Line 8283  void prevbackforecast(char fileres[], do Line 8326  void prevbackforecast(char fileres[], do
   if(estepm < stepm){    if(estepm < stepm){
     printf ("Problem %d lower than %d\n",estepm, stepm);      printf ("Problem %d lower than %d\n",estepm, stepm);
   }    }
   else  hstepm=estepm;    else{
       hstepm=estepm;   
     }
     if(estepm >= stepm){ /* Yes every two year */
       stepsize=2;
     }
       
   hstepm=hstepm/stepm;    hstepm=hstepm/stepm;
   yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and    yp1=modf(dateintmean,&yp);/* extracts integral of datemean in yp  and
Line 8304  void prevbackforecast(char fileres[], do Line 8352  void prevbackforecast(char fileres[], do
       
   fprintf(ficresfb,"#****** Routine prevbackforecast **\n");    fprintf(ficresfb,"#****** Routine prevbackforecast **\n");
       
   /*          if (h==(int)(YEARM*yearp)){ */  
   /* for(cptcov=1, k=0;cptcov<=i1;cptcov++){ */  
   /*   for(cptcod=1;cptcod<=ncodemax[cptcoveff];cptcod++){ */  
   /*     k=k+1; */  
   for(nres=1; nres <= nresult; nres++) /* For each resultline */    for(nres=1; nres <= nresult; nres++) /* For each resultline */
   for(k=1; k<=i1;k++){    for(k=1; k<=i1;k++){
     if(i1 != 1 && TKresult[nres]!= k)      if(i1 != 1 && TKresult[nres]!= k)
Line 8333  void prevbackforecast(char fileres[], do Line 8377  void prevbackforecast(char fileres[], do
       /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {  */        /* for (yearp=0; yearp<=(anproj2-anproj1);yearp +=stepsize) {  */
       fprintf(ficresfb,"\n");        fprintf(ficresfb,"\n");
       fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);        fprintf(ficresfb,"\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);
       printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp);        /* printf("\n# Back Forecasting at date %.lf/%.lf/%.lf ",jback1,mback1,anback1+yearp); */
       /* for (agec=fage; agec>=(ageminpar-1); agec--){ */        /* for (agec=bage; agec<=agemax-1; agec++){  /\* testing *\/ */
       /*        nhstepm=(int) rint((agelim-agec)*YEARM/stepm); */        for (agec=bage; agec<=fage; agec++){  /* testing */
       for (agec=bage; agec<=agemax-1; agec++){  /* testing */  
         /* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/          /* We compute bij at age agec over nhstepm, nhstepm decreases when agec increases because of agemax;*/
         nhstepm=(int) rint((agec-agelim)*YEARM/stepm);          nhstepm=(int) (agec-agelim) *YEARM/stepm;/*     nhstepm=(int) rint((agec-agelim)*YEARM/stepm);*/
         nhstepm = nhstepm/hstepm;          nhstepm = nhstepm/hstepm;
         p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);          p3mat=ma3x(1,nlstate+ndeath,1, nlstate+ndeath, 0,nhstepm);
         oldm=oldms;savm=savms;          oldm=oldms;savm=savms;
         /* computes hbxij at age agec over 1 to nhstepm */          /* computes hbxij at age agec over 1 to nhstepm */
           /* printf("####prevbackforecast debug  agec=%.2f nhstepm=%d\n",agec, nhstepm);fflush(stdout); */
         hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);          hbxij(p3mat,nhstepm,agec,hstepm,p,prevacurrent,nlstate,stepm, k, nres);
         /* hpxij(p3mat,nhstepm,agec,hstepm,p,             nlstate,stepm,oldm,savm, k,nres); */          /* hpxij(p3mat,nhstepm,agec,hstepm,p,             nlstate,stepm,oldm,savm, k,nres); */
         /* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */          /* Then we print p3mat for h corresponding to the right agec+h*stepms=yearp */
Line 8360  void prevbackforecast(char fileres[], do Line 8404  void prevbackforecast(char fileres[], do
           ppij=0.;ppi=0.;            ppij=0.;ppi=0.;
           for(j=1; j<=nlstate;j++) {            for(j=1; j<=nlstate;j++) {
             /* if (mobilav==1) */              /* if (mobilav==1) */
             ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k];              ppij=ppij+p3mat[i][j][h]*prevacurrent[(int)agec][j][k];
             ppi=ppi+mobaverage[(int)agec][j][k];              ppi=ppi+prevacurrent[(int)agec][j][k];
               /* ppij=ppij+p3mat[i][j][h]*mobaverage[(int)agec][j][k]; */
               /* ppi=ppi+mobaverage[(int)agec][j][k]; */
               /* else { */                /* else { */
               /*        ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */                /*        ppij=ppij+p3mat[i][j][h]*probs[(int)(agec)][i][k]; */
               /* } */                /* } */
Line 8386  void prevbackforecast(char fileres[], do Line 8432  void prevbackforecast(char fileres[], do
                   
 }  }
   
   /* Variance of prevalence limit: varprlim */
    void varprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **prlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){
       /*------- Variance of period (stable) prevalence------*/   
    
      char fileresvpl[FILENAMELENGTH];  
      FILE *ficresvpl;
      double **oldm, **savm;
      double **varpl; /* Variances of prevalence limits by age */   
      int i1, k, nres, j ;
      
       strcpy(fileresvpl,"VPL_");
       strcat(fileresvpl,fileresu);
       if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {
         printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);
         exit(0);
       }
       printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);
       fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);
       
       /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){
         for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/
       
       i1=pow(2,cptcoveff);
       if (cptcovn < 1){i1=1;}
   
       for(nres=1; nres <= nresult; nres++) /* For each resultline */
       for(k=1; k<=i1;k++){
         if(i1 != 1 && TKresult[nres]!= k)
           continue;
         fprintf(ficresvpl,"\n#****** ");
         printf("\n#****** ");
         fprintf(ficlog,"\n#****** ");
         for(j=1;j<=cptcoveff;j++) {
           fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
           fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
           printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
         }
         for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
           printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
           fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
           fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
         } 
         fprintf(ficresvpl,"******\n");
         printf("******\n");
         fprintf(ficlog,"******\n");
         
         varpl=matrix(1,nlstate,(int) bage, (int) fage);
         oldm=oldms;savm=savms;
         varprevlim(fileresvpl, ficresvpl, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, ncvyearp, k, strstart, nres);
         free_matrix(varpl,1,nlstate,(int) bage, (int)fage);
         /*}*/
       }
       
       fclose(ficresvpl);
       printf("done variance-covariance of period prevalence\n");fflush(stdout);
       fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);
   
    }
   /* Variance of back prevalence: varbprlim */
    void varbprlim(char fileresu[], int nresult, double ***prevacurrent, int mobilavproj, double bage, double fage, double **bprlim, int *ncvyearp, double ftolpl, double p[], double **matcov, double *delti, int stepm, int cptcoveff){
         /*------- Variance of back (stable) prevalence------*/
   
      char fileresvbl[FILENAMELENGTH];  
      FILE  *ficresvbl;
   
      double **oldm, **savm;
      double **varbpl; /* Variances of back prevalence limits by age */   
      int i1, k, nres, j ;
   
      strcpy(fileresvbl,"VBL_");
      strcat(fileresvbl,fileresu);
      if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {
        printf("Problem with variance of back (stable) prevalence  resultfile: %s\n", fileresvbl);
        exit(0);
      }
      printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);
      fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);
      
      
      i1=pow(2,cptcoveff);
      if (cptcovn < 1){i1=1;}
      
      for(nres=1; nres <= nresult; nres++) /* For each resultline */
        for(k=1; k<=i1;k++){
          if(i1 != 1 && TKresult[nres]!= k)
            continue;
          fprintf(ficresvbl,"\n#****** ");
          printf("\n#****** ");
          fprintf(ficlog,"\n#****** ");
          for(j=1;j<=cptcoveff;j++) {
            fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
            fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
            printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);
          }
          for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */
            printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
            fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
            fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);
          }
          fprintf(ficresvbl,"******\n");
          printf("******\n");
          fprintf(ficlog,"******\n");
          
          varbpl=matrix(1,nlstate,(int) bage, (int) fage);
          oldm=oldms;savm=savms;
          
          varbrevlim(fileresvbl, ficresvbl, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, ncvyearp, k, strstart, nres);
          free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);
          /*}*/
        }
      
      fclose(ficresvbl);
      printf("done variance-covariance of back prevalence\n");fflush(stdout);
      fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);
   
    } /* End of varbprlim */
   
 /************** Forecasting *****not tested NB*************/  /************** Forecasting *****not tested NB*************/
 /* void populforecast(char fileres[], double anpyram,double mpyram,double jpyram,double ageminpar, double agemax,double dateprev1, double dateprev2s, 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 dateprev2s, int mobilav, double agedeb, double fage, int popforecast, char popfile[], double anpyram1,double p[], int i2){ */
       
Line 10443  int main(int argc, char *argv[]) Line 10606  int main(int argc, char *argv[])
   int nres=0;    int nres=0;
   int endishere=0;    int endishere=0;
   
     int ncurrv=0; /* Temporary variable */
     
   char ca[32], cb[32];    char ca[32], cb[32];
   /*  FILE *fichtm; *//* Html File */    /*  FILE *fichtm; *//* Html File */
   /* FILE *ficgp;*/ /*Gnuplot File */    /* FILE *ficgp;*/ /*Gnuplot File */
Line 10494  int main(int argc, char *argv[]) Line 10659  int main(int argc, char *argv[])
   double *delti; /* Scale */    double *delti; /* Scale */
   double ***eij, ***vareij;    double ***eij, ***vareij;
   double **varpl; /* Variances of prevalence limits by age */    double **varpl; /* Variances of prevalence limits by age */
   double **varbpl; /* Variances of back prevalence limits by age */  
   double *epj, vepp;    double *epj, vepp;
   
   double dateprev1, dateprev2,jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000;    double dateprev1, dateprev2;
   double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000;    double jproj1=1,mproj1=1,anproj1=2000,jproj2=1,mproj2=1,anproj2=2000, dateproj1=0, dateproj2=0;
     double jback1=1,mback1=1,anback1=2000,jback2=1,mback2=1,anback2=2000, dateback1=0, dateback2=0;
   
   double **ximort;    double **ximort;
   char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";    char *alph[]={"a","a","b","c","d","e"}, str[4]="1234";
Line 11326  Title=%s <br>Datafile=%s Firstpass=%d La Line 11492  Title=%s <br>Datafile=%s Firstpass=%d La
               firstpass, lastpass,  stepm,  weightopt, model);                firstpass, lastpass,  stepm,  weightopt, model);
   
   fprintf(fichtm,"\n");    fprintf(fichtm,"\n");
   fprintf(fichtm,"<br>Total number of observations=%d <br>\n\    fprintf(fichtm,"<h4>Parameter line 2</h4><ul><li>Tolerance for the convergence of the likelihood: ftol=%f \n<li>Interval for the elementary matrix (in month): stepm=%d",\
             ftol, stepm);
     fprintf(fichtm,"\n<li>Number of fixed dummy covariates: ncovcol=%d ", ncovcol);
     ncurrv=1;
     for(i=ncurrv; i <=ncovcol; i++) fprintf(fichtm,"V%d ", i);
     fprintf(fichtm,"\n<li> Number of fixed quantitative variables: nqv=%d ", nqv); 
     ncurrv=i;
     for(i=ncurrv; i <=ncurrv-1+nqv; i++) fprintf(fichtm,"V%d ", i);
     fprintf(fichtm,"\n<li> Number of time varying (wave varying) covariates: ntv=%d ", ntv);
     ncurrv=i;
     for(i=ncurrv; i <=ncurrv-1+ntv; i++) fprintf(fichtm,"V%d ", i);
     fprintf(fichtm,"\n<li>Number of quantitative time varying covariates: nqtv=%d ", nqtv);
     ncurrv=i;
     for(i=ncurrv; i <=ncurrv-1+nqtv; i++) fprintf(fichtm,"V%d ", i);
     fprintf(fichtm,"\n<li>Weights column \n<br>Number of alive states: nlstate=%d <br>Number of death states (not really implemented): ndeath=%d \n<li>Number of waves: maxwav=%d \n<li>Parameter for maximization (1), using parameter values (0), for design of parameters and variance-covariance matrix: mle=%d \n<li>Does the weight column be taken into account (1), or not (0): weight=%d</ul>\n", \
              nlstate, ndeath, maxwav, mle, weightopt);
   
     fprintf(fichtm,"<h4> Diagram of states <a href=\"%s_.svg\">%s_.svg</a></h4> \n\
   <img src=\"%s_.svg\">", subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"),subdirf2(optionfilefiname,"D_"));
   
     
     fprintf(fichtm,"\n<h4>Some descriptive statistics </h4>\n<br>Total number of observations=%d <br>\n\
 Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\  Youngest age at first (selected) pass %.2f, oldest age %.2f<br>\n\
 Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\  Interval (in months) between two waves: Min=%d Max=%d Mean=%.2lf<br>\n",\
           imx,agemin,agemax,jmin,jmax,jmean);    imx,agemin,agemax,jmin,jmax,jmean);
   pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
   oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
   newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */    newms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
Line 11889  Please run with mle=-1 to get a correct Line 12076  Please run with mle=-1 to get a correct
           fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);            fprintf(ficlog,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);
           fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,mobilavproj);            fprintf(ficres,"prevforecast=%d starting-proj-date=%.lf/%.lf/%.lf final-proj-date=%.lf/%.lf/%.lf mobil_average=%d\n",prevfcast,jproj1,mproj1,anproj1,jproj2,mproj2,anproj2,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.*/
             dateproj1=anproj1+(mproj1-1)/12.+(jproj1-1)/365.;
             dateproj2=anproj2+(mproj2-1)/12.+(jproj2-1)/365.;
   
         }          }
         break;          break;
       case 12:        case 12:
Line 11904  Please run with mle=-1 to get a correct Line 12094  Please run with mle=-1 to get a correct
           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.*/
             dateback1=anback1+(mback1-1)/12.+(jback1-1)/365.;
             dateback2=anback2+(mback2-1)/12.+(jback2-1)/365.;
         }          }
         break;          break;
       case 13:        case 13:
Line 11954  Please run with mle=-1 to get a correct Line 12146  Please run with mle=-1 to get a correct
 This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\  This is probably because your parameter file doesn't \n  contain the exact number of lines (or columns) corresponding to your model line.\n\
 Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);  Please run with mle=-1 to get a correct covariance matrix.\n",ageminpar,agemaxpar);
     }else{      }else{
       printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1);        /* printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,fage, prevfcast, backcast, pathc,p, (int)anproj1-(int)agemin, (int)anback1-(int)agemax+1); */
         printinggnuplot(fileresu, optionfilefiname,ageminpar,agemaxpar,bage, fage, prevfcast, backcast, pathc,p, (int)anproj1-bage, (int)anback1-fage);
     }      }
     printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \      printinghtml(fileresu,title,datafile, firstpass, lastpass, stepm, weightopt, \
                  model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \                   model,imx,jmin,jmax,jmean,rfileres,popforecast,mobilav,prevfcast,mobilavproj,backcast, estepm, \
                  jprev1,mprev1,anprev1,dateprev1,jprev2,mprev2,anprev2,dateprev2);                   jprev1,mprev1,anprev1,dateprev1, dateproj1, dateback1,jprev2,mprev2,anprev2,dateprev2,dateproj2, dateback2);
                                   
     /*------------ free_vector  -------------*/      /*------------ free_vector  -------------*/
     /*  chdir(path); */      /*  chdir(path); */
Line 11994  Please run with mle=-1 to get a correct Line 12187  Please run with mle=-1 to get a correct
     k=1;      k=1;
     varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);      varprob(optionfilefiname, matcov, p, delti, nlstate, bage, fage,k,Tvar,nbcode, ncodemax,strstart);
           
     /* Prevalence for each covariates in probs[age][status][cov] */      /* Prevalence for each covariate combination in probs[age][status][cov] */
     probs= ma3x(1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);      probs= ma3x(AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
     for(i=1;i<=AGESUP;i++)      for(i=AGEINF;i<=AGESUP;i++)
       for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */        for(j=1;j<=nlstate+ndeath;j++) /* ndeath is useless but a necessity to be compared with mobaverages */
         for(k=1;k<=ncovcombmax;k++)          for(k=1;k<=ncovcombmax;k++)
           probs[i][j][k]=0.;            probs[i][j][k]=0.;
     prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);      prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, 
                  ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass);
     if (mobilav!=0 ||mobilavproj !=0 ) {      if (mobilav!=0 ||mobilavproj !=0 ) {
       mobaverages= ma3x(1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);        mobaverages= ma3x(AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
       for(i=1;i<=AGESUP;i++)        for(i=AGEINF;i<=AGESUP;i++)
         for(j=1;j<=nlstate+ndeath;j++)          for(j=1;j<=nlstate+ndeath;j++)
           for(k=1;k<=ncovcombmax;k++)            for(k=1;k<=ncovcombmax;k++)
             mobaverages[i][j][k]=0.;              mobaverages[i][j][k]=0.;
Line 12015  Please run with mle=-1 to get a correct Line 12209  Please run with mle=-1 to get a correct
           fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);            fprintf(ficlog," Error in movingaverage mobilav=%d\n",mobilav);
           printf(" Error in movingaverage mobilav=%d\n",mobilav);            printf(" Error in movingaverage mobilav=%d\n",mobilav);
         }          }
       }        } else if (mobilavproj !=0) {
       /* else if(mobilavproj==-1){ /\* Forcing raw observed prevalences *\/ */  
       /*        for(i=1;i<=AGESUP;i++) */  
       /*          for(j=1;j<=nlstate;j++) */  
       /*            for(k=1;k<=ncovcombmax;k++) */  
       /*              mobaverages[i][j][k]=probs[i][j][k]; */  
       /*        /\* /\\* Prevalence for each covariates in probs[age][status][cov] *\\/ *\/ */  
       /*        /\* prevalence(probs, ageminpar, agemaxpar, s, agev, nlstate, imx, Tvar, nbcode, ncodemax, mint, anint, dateprev1, dateprev2, firstpass, lastpass); *\/ */  
       /* } */  
       else if (mobilavproj !=0) {  
         printf("Movingaveraging projected observed prevalence\n");          printf("Movingaveraging projected observed prevalence\n");
         fprintf(ficlog,"Movingaveraging projected observed prevalence\n");          fprintf(ficlog,"Movingaveraging projected observed prevalence\n");
         if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){          if (movingaverage(probs, ageminpar, agemaxpar, mobaverage, mobilavproj)!=0){
           fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);            fprintf(ficlog," Error in movingaverage mobilavproj=%d\n",mobilavproj);
           printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);            printf(" Error in movingaverage mobilavproj=%d\n",mobilavproj);
         }          }
         }else{
           printf("Internal error moving average\n");
           fflush(stdout);
           exit(1);
       }        }
     }/* end if moving average */      }/* end if moving average */
           
     /*---------- Forecasting ------------------*/      /*---------- Forecasting ------------------*/
     /*if((stepm == 1) && (strcmp(model,".")==0)){*/  
     if(prevfcast==1){      if(prevfcast==1){
       /*    if(stepm ==1){*/        /*    if(stepm ==1){*/
       prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);        prevforecast(fileresu, anproj1, mproj1, jproj1, agemin, agemax, dateprev1, dateprev2, mobilavproj, mobaverage, bage, fage, firstpass, lastpass, anproj2, p, cptcoveff);
     }      }
   
       /* Backcasting */
     if(backcast==1){      if(backcast==1){
       ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);                ddnewms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);        
       ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);                ddoldms=matrix(1,nlstate+ndeath,1,nlstate+ndeath);        
Line 12048  Please run with mle=-1 to get a correct Line 12238  Please run with mle=-1 to get a correct
       /*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/        /*--------------- Back Prevalence limit  (period or stable prevalence) --------------*/
   
       bprlim=matrix(1,nlstate,1,nlstate);        bprlim=matrix(1,nlstate,1,nlstate);
   
       back_prevalence_limit(p, bprlim,  ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);        back_prevalence_limit(p, bprlim,  ageminpar, agemaxpar, ftolpl, &ncvyear, dateprev1, dateprev2, firstpass, lastpass, mobilavproj);
       fclose(ficresplb);        fclose(ficresplb);
   
       hBijx(p, bage, fage, mobaverage);        hBijx(p, bage, fage, mobaverage);
       fclose(ficrespijb);        fclose(ficrespijb);
       /* free_matrix(bprlim,1,nlstate,1,nlstate); /\*here or after loop ? *\/ */  
   
       prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2, mobilavproj,        prevbackforecast(fileresu, mobaverage, anback1, mback1, jback1, agemin, agemax, dateprev1, dateprev2,
          bage, fage, firstpass, lastpass, anback2, p, cptcoveff);                         mobilavproj, bage, fage, firstpass, lastpass, anback2, p, cptcoveff);
         varbprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, bprlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
   
       /*------- Variance of back (stable) prevalence------*/     
               
       strcpy(fileresvbl,"VBL_");        free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */
       strcat(fileresvbl,fileresu);  
       if((ficresvbl=fopen(fileresvbl,"w"))==NULL) {  
         printf("Problem with variance of back (stable) prevalence  resultfile: %s\n", fileresvbl);  
         exit(0);  
       }  
       printf("Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(stdout);  
       fprintf(ficlog, "Computing Variance-covariance of back (stable) prevalence: file '%s' ...", fileresvbl);fflush(ficlog);  
         
       /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){  
         for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/  
         
       i1=pow(2,cptcoveff);  
       if (cptcovn < 1){i1=1;}  
         
       for(nres=1; nres <= nresult; nres++) /* For each resultline */  
         for(k=1; k<=i1;k++){  
           if(i1 != 1 && TKresult[nres]!= k)  
             continue;  
           fprintf(ficresvbl,"\n#****** ");  
           printf("\n#****** ");  
           fprintf(ficlog,"\n#****** ");  
           for(j=1;j<=cptcoveff;j++) {  
             fprintf(ficresvbl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);  
             fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);  
             printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);  
           }  
           for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */  
             printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);  
             fprintf(ficresvbl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);  
             fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);  
           }       
           fprintf(ficresvbl,"******\n");  
           printf("******\n");  
           fprintf(ficlog,"******\n");  
             
           varbpl=matrix(1,nlstate,(int) bage, (int) fage);  
           oldm=oldms;savm=savms;  
             
           varbrevlim(fileres, varbpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, bprlim, ftolpl, mobilavproj, &ncvyear, k, strstart, nres);  
           free_matrix(varbpl,1,nlstate,(int) bage, (int)fage);  
           /*}*/  
         }  
       
       fclose(ficresvbl);  
       printf("done variance-covariance of back prevalence\n");fflush(stdout);  
       fprintf(ficlog,"done variance-covariance of back prevalence\n");fflush(ficlog);  
   
       free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);        free_matrix(ddnewms, 1, nlstate+ndeath, 1, nlstate+ndeath);
       free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);        free_matrix(ddsavms, 1, nlstate+ndeath, 1, nlstate+ndeath);
       free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);        free_matrix(ddoldms, 1, nlstate+ndeath, 1, nlstate+ndeath);
     }      }    /* end  Backcasting */
     free_matrix(bprlim,1,nlstate,1,nlstate); /*here or after loop ? */  
     
     
     /* ------ Other prevalence ratios------------ */      /* ------ Other prevalence ratios------------ */
Line 12165  Please run with mle=-1 to get a correct Line 12307  Please run with mle=-1 to get a correct
     fclose(ficreseij);      fclose(ficreseij);
     printf("done evsij\n");fflush(stdout);      printf("done evsij\n");fflush(stdout);
     fprintf(ficlog,"done evsij\n");fflush(ficlog);      fprintf(ficlog,"done evsij\n");fflush(ficlog);
   
                                   
     /*---------- State-specific expectancies and variances ------------*/      /*---------- State-specific expectancies and variances ------------*/
                                   
                   
     strcpy(filerest,"T_");      strcpy(filerest,"T_");
     strcat(filerest,fileresu);      strcat(filerest,fileresu);
     if((ficrest=fopen(filerest,"w"))==NULL) {      if((ficrest=fopen(filerest,"w"))==NULL) {
Line 12177  Please run with mle=-1 to get a correct Line 12319  Please run with mle=-1 to get a correct
     }      }
     printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);      printf("Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(stdout);
     fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);      fprintf(ficlog,"Computing Total Life expectancies with their standard errors: file '%s' ...\n", filerest); fflush(ficlog);
                   
   
     strcpy(fileresstde,"STDE_");      strcpy(fileresstde,"STDE_");
     strcat(fileresstde,fileresu);      strcat(fileresstde,fileresu);
     if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {      if((ficresstdeij=fopen(fileresstde,"w"))==NULL) {
Line 12206  Please run with mle=-1 to get a correct Line 12346  Please run with mle=-1 to get a correct
     printf("      Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout);      printf("      Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(stdout);
     fprintf(ficlog,"      Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog);      fprintf(ficlog,"      Computing Variance-covariance of State-specific Expectancies: file '%s' ... ", fileresv);fflush(ficlog);
   
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){  
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/  
             
     i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */      i1=pow(2,cptcoveff); /* Number of combination of dummy covariates */
     if (cptcovn < 1){i1=1;}      if (cptcovn < 1){i1=1;}
           
Line 12270  Please run with mle=-1 to get a correct Line 12407  Please run with mle=-1 to get a correct
       vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);        vareij=ma3x(1,nlstate,1,nlstate,(int) bage, (int) fage);
       pstamp(ficrest);        pstamp(ficrest);
               
               epj=vector(1,nlstate+1);
       for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/        for(vpopbased=0; vpopbased <= popbased; vpopbased++){ /* Done for vpopbased=0 and vpopbased=1 if popbased==1*/
         oldm=oldms;savm=savms; /* ZZ Segmentation fault */          oldm=oldms;savm=savms; /* ZZ Segmentation fault */
         cptcod= 0; /* To be deleted */          cptcod= 0; /* To be deleted */
Line 12286  Please run with mle=-1 to get a correct Line 12423  Please run with mle=-1 to get a correct
         for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);          for (i=1;i<=nlstate;i++) fprintf(ficrest,"e.%d (std) ",i);
         fprintf(ficrest,"\n");          fprintf(ficrest,"\n");
         /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */          /* printf("Which p?\n"); for(i=1;i<=npar;i++)printf("p[i=%d]=%lf,",i,p[i]);printf("\n"); */
         epj=vector(1,nlstate+1);  
         printf("Computing age specific period (stable) prevalences in each health state \n");          printf("Computing age specific period (stable) prevalences in each health state \n");
         fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");          fprintf(ficlog,"Computing age specific period (stable) prevalences in each health state \n");
         for(age=bage; age <=fage ;age++){          for(age=bage; age <=fage ;age++){
Line 12324  Please run with mle=-1 to get a correct Line 12460  Please run with mle=-1 to get a correct
           fprintf(ficrest,"\n");            fprintf(ficrest,"\n");
         }          }
       } /* End vpopbased */        } /* End vpopbased */
         free_vector(epj,1,nlstate+1);
       free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);        free_ma3x(eij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);        free_ma3x(vareij,1,nlstate,1,nlstate,(int) bage, (int)fage);
       free_vector(epj,1,nlstate+1);  
       printf("done selection\n");fflush(stdout);        printf("done selection\n");fflush(stdout);
       fprintf(ficlog,"done selection\n");fflush(ficlog);        fprintf(ficlog,"done selection\n");fflush(ficlog);
               
       /*}*/  
     } /* End k selection */      } /* End k selection */
   
     printf("done State-specific expectancies\n");fflush(stdout);      printf("done State-specific expectancies\n");fflush(stdout);
     fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);      fprintf(ficlog,"done State-specific expectancies\n");fflush(ficlog);
   
     /*------- Variance of period (stable) prevalence------*/         /* variance-covariance of period prevalence*/
           varprlim(fileresu, nresult, mobaverage, mobilavproj, bage, fage, prlim, &ncvyear, ftolpl, p, matcov, delti, stepm, cptcoveff);
     strcpy(fileresvpl,"VPL_");  
     strcat(fileresvpl,fileresu);  
     if((ficresvpl=fopen(fileresvpl,"w"))==NULL) {  
       printf("Problem with variance of period (stable) prevalence  resultfile: %s\n", fileresvpl);  
       exit(0);  
     }  
     printf("Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(stdout);  
     fprintf(ficlog, "Computing Variance-covariance of period (stable) prevalence: file '%s' ...", fileresvpl);fflush(ficlog);  
       
     /*for(cptcov=1,k=0;cptcov<=i1;cptcov++){  
       for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){*/  
       
     i1=pow(2,cptcoveff);  
     if (cptcovn < 1){i1=1;}  
   
     for(nres=1; nres <= nresult; nres++) /* For each resultline */  
     for(k=1; k<=i1;k++){  
       if(i1 != 1 && TKresult[nres]!= k)  
         continue;  
       fprintf(ficresvpl,"\n#****** ");  
       printf("\n#****** ");  
       fprintf(ficlog,"\n#****** ");  
       for(j=1;j<=cptcoveff;j++) {  
         fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);  
         fprintf(ficlog,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);  
         printf("V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtabm(k,j)]);  
       }  
       for (j=1; j<= nsq; j++){ /* For each selected (single) quantitative value */  
         printf(" V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);  
         fprintf(ficresvpl," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);  
         fprintf(ficlog," V%d=%f ",Tvqresult[nres][j],Tqresult[nres][j]);  
       }   
       fprintf(ficresvpl,"******\n");  
       printf("******\n");  
       fprintf(ficlog,"******\n");  
         
       varpl=matrix(1,nlstate,(int) bage, (int) fage);  
       oldm=oldms;savm=savms;  
       varprevlim(fileres, varpl, matcov, p, delti, nlstate, stepm, (int) bage, (int) fage, oldm, savm, prlim, ftolpl, &ncvyear, k, strstart, nres);  
       free_matrix(varpl,1,nlstate,(int) bage, (int)fage);  
       /*}*/  
     }  
       
     fclose(ficresvpl);  
     printf("done variance-covariance of period prevalence\n");fflush(stdout);  
     fprintf(ficlog,"done variance-covariance of period prevalence\n");fflush(ficlog);  
   
           
     free_vector(weight,1,n);      free_vector(weight,1,n);
Line 12402  Please run with mle=-1 to get a correct Line 12491  Please run with mle=-1 to get a correct
           
     /*---------- End : free ----------------*/      /*---------- End : free ----------------*/
     if (mobilav!=0 ||mobilavproj !=0)      if (mobilav!=0 ||mobilavproj !=0)
       free_ma3x(mobaverages,1, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */        free_ma3x(mobaverages,AGEINF, AGESUP,1,nlstate+ndeath, 1,ncovcombmax); /* We need to have a squared matrix with prevalence of the dead! */
     free_ma3x(probs,1,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);      free_ma3x(probs,AGEINF,AGESUP,1,nlstate+ndeath, 1,ncovcombmax);
     free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */      free_matrix(prlim,1,nlstate,1,nlstate); /*here or after loop ? */
     free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);      free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
   }  /* mle==-3 arrives here for freeing */    }  /* mle==-3 arrives here for freeing */
Line 12420  Please run with mle=-1 to get a correct Line 12509  Please run with mle=-1 to get a correct
   /*free_vector(delti,1,npar);*/    /*free_vector(delti,1,npar);*/
   free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);     free_ma3x(delti3,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel); 
   free_matrix(agev,1,maxwav,1,imx);    free_matrix(agev,1,maxwav,1,imx);
     free_ma3x(paramstart,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
   free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);    free_ma3x(param,1,nlstate,1, nlstate+ndeath-1,1,ncovmodel);
       
   free_ivector(ncodemax,1,NCOVMAX);    free_ivector(ncodemax,1,NCOVMAX);
Removed from v.1.268  
changed lines
  Added in v.1.275

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